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ABSTRACT 


Frequently,  in  structural  system  identification  (model  updating  or  damage  detection),  the 
available  set  of  data  is  incomplete,  both  spatially  and  in  modal  content.  This 
incompleteness  leads  to  the  solution  of  underdetennined  linear  systems.  In  order  to 
improve  the  identification,  additional  independent  measured  data  must  be  found.  In  the 
past,  it  has  been  shown  that  such  data  can  be  easily  obtained  from  the  application  of 
Artificial  Boundary  Conditions  (ABC),  imposed  on  both  the  baseline  FE  models  and  the 
measured  frequency  response  data.  This  can  be  accomplished  without  any  physical 
modifications  to  the  experiment  and,  hence,  no  additional  expense  on  different  systems, 
or  more  than  once,  in  order  to  get  the  modal  data  needed  for  the  analysis.  In  this  thesis, 
the  procedure  of  sensitivity-based  structural  system  identification,  using  ABCs,  and 
enhanced  by  parameter  grouping/clustering,  is  examined.  It  is  shown  that  the  optimal 
sensitivity  matrix  is  a  square  and  diagonal  dominant  one,  which  can  be  used  with  quite 
accurate  results  both  for  localization  of  parameter  errors,  and  the  determination  of  the 
magnitude  of  parameter  error.  The  numerous  ABC  configurations  available,  even  from  a 
small  measured  data  set,  allow  an  optimal  sensitivity  matrix  to  be  found  for  many 
parameters.  These  concepts  are  demonstrated  using  simulated  measurements  along  with 
finite  element  models. 
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I.  INTRODUCTION 


A.  MODEL  UPDATING  AND  DAMAGE  DETECTION 

Finite  element  analysis  (FE)  is  a  computational  tool  that  is  widely  used  in  cases 
where  a  closed-form  solution  to  a  partial  differential  equation  is  not  available.  In 
structural  dynamics,  FE  models  are  used  to  predict  the  dynamic  response  of  a  system  to 
external  stimuli,  such  as  dynamic  excitations,  or  to  assess  the  effect  of  parameter 
changes.  The  term  “parameter”  stands  for  all  the  individual  physical  and  geometric  data 
that  define  the  structure,  such  as  geometry,  elastic  moduli  and  element  masses/densities. 

The  use  of  FE  models  aims  at  achieving  higher  accuracy  in  terms  of  response 
prediction;  however,  doing  so  may  not  always  be  possible  for  several  reasons.  Some  of 
them  are  usually  connected  to  uncertainties  regarding  the  governing  equation  of  the 
system,  to  assumptions  made  during  the  derivation,  to  inaccurate  physical  boundary 
conditions,  to  inexact  characterization  of  joints  and  connections,  etc.  All  these 
discrepancies  are  revealed  and  attempted  to  be  overcome,  when  a  necessary  procedure  for 
the  validation  of  the  FE  model  takes  place.  This  procedure  is  called  FEM  updating.  It 
involves  measuring  structure  response,  identifying  modal  data,  comparing  the  modal 
parameters  to  those  extracted  from  the  FE  model  and  finally,  modifying  appropriately  the 
FE  model  to  match  the  system’s  design  requirements.  FEM  updating  is  of  vital 
importance  as  it  verifies  that  the  model  can  be  used  with  acceptable  reliability  for  a 
number  of  research  fields  in  modal  analysis,  such  as  damage  detection  and  error 
localization. 

B.  BACKGROUND 

Over  the  last  decades,  a  large  volume  of  literature  on  structural  identification, 
including  both  model  updating  and  error  localization  has  been  produced.  Friswell  and 
Mottershead  in  their  survey  [1],  extensively  discuss  methods  that  deal  with  various 
aspects  of  this  field.  Most  of  them  are  eigenvalue  sensitivity-based  and  are  thoroughly 
presented  and  analyzed  in  [2].  Each  developed  method  has,  as  a  main  objective,  to 
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overcome  problems  that  arise  during  the  procedure,  which  in  general,  converges  to  the 
solving  of  a  linear  system.  The  most  common  issue  with  the  solution  is  the  ill- 
conditioning  of  these  systems.  This  ill-condition  results  from  the  lack  of  input  data, 
which  frequently  leads  to  a  severely  underdetennined  system  of  equations.  The  lack  of 
input  data  is  mainly  caused  either  by  the  geometric  limitation  of  measurement  locations 
on  a  structure  (spatial  incompleteness)  or  by  the  bounds  of  the  available  bandwidth 
during  modal  testing  (modal  incompleteness).  It  is  often  observed  that  ill-conditioning 
may  be  a  concern  due  to  dependencies  of  the  sensitivity  matrix,  or  even  due  to  the 
physical  geometry  of  the  structure,  in  cases  where  symmetry  is  present. 

Dealing  with  this  problem,  Lallement  and  Piranda  in  [3]  are  the  first  to  introduce 
the  development  of  a  rational  procedure  for  parameter  selection — which  appears  to  be  the 
main  core  in  many  updating  methods — to  facilitate  the  improvement  of  the  conditioning 
of  the  equations.  This  has  to  do  with  an  iterative  procedure  to  find  an  appropriate  subset 
of  parameters  that  produce  an  optimal  solution  commonly  known  as  forward  selection. 
Then,  the  importance  of  vector  subspaces  and  their  wide  application  range  in  modal 
analysis  is  revealed  [4],  and  in  [5]  the  concept  of  selecting  parameters  groups  using 
angles  between  subspaces  is  implemented,  also  incorporating  side  constraints  to  the 
equations  to  be  solved.  The  methods  of  subset  selection  are  outlined  and  extended  to  the 
selection  of  groups  of  parameters.  Reference  [6]  is  a  representative  example  of  the 
selection  and  updating  of  parameters  in  an  actual  application.  Another  approach  to  the 
same  objective,  parameter  selection,  is  offered  in  [7]  where  an  automated  procedure  to 
guide  the  parameter  selection  is  suggested  based  on  simple  observations.  The  idea  of 
grouping  element  sensitivities  together,  under  the  condition  of  selecting  as  few 
parameters  as  possible  that  comply  with  the  above  stated  observations,  is  presented. 
These  observations  deal  with  the  fact  that  the  number  of  updating  parameters  can  be 
reduced,  but  only  by  sacrificing  the  total  sensitivity,  which  is  a  concept  that  will  be 
covered  and  used  in  Chapter  IV  of  this  thesis.  References  [8],  [9]  and  [10]  address  this 
concept  in  a  different  manner,  using  the  principles  of  cluster  analysis,  along  with  optimal 
subset  selection,  and  show  how  these  methods  can  be  successfully  used  in  a  number  of 
applications. 
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Another  path  towards  the  improving  the  ill-conditioning  aims  at  the  enrichment  of 
the  available  data,  in  order  to  decrease  incompleteness,  or  be  less  underdetennined.  Of 
course,  care  has  to  be  exercised  in  preventing  the  addition  of  redundant  data  and  in 
accordance  with  feasibility  conditions.  A  first  step  of  doing  so  without  requiring  any 
modifications  to  the  physical  structure  is  made  by  Gordis  [11],  where  the  concept  of 
artificial  constraints  is  initially  introduced.  The  term  artificial  refers  to  the  fact  that  a 
constraint  is  applied  computationally  to  the  measured  data,  as  well  as  to  the  FE  model. 
Another  similar  reference  is  made  by  Rade,  Lallement  and  Silva  in  [12],  stated  here  as 
Modified  Configurations,  describing  related  ways  different  artificial  systems  are  used  to 
enlarge  the  experimental  data.  In  this  work,  the  criterion  established  is  the  minimization 
of  the  condition  number  of  the  sensitivity  matrix,  or  in  other  words,  the  maximization  of 
the  degree  of  linear  independence  between  the  columns  and  the  rows  of  the  sensitivity 
matrix.  A  theoretical  investigation  of  this  concept  is  made  in  [13],  where  the  term  ABC, 
standing  for  Artificial  Boundary  Condition,  is  first  introduced.  It  proves  that,  with  no 
additional  testing,  all  the  natural  frequencies  corresponding  to  ABC  systems  and  the 
needed  additional  data  can  be  obtained  from  any  square  FRF  matrix  measured  in  a  modal 
test.  In  addition,  it  shows  that  structural  sensitivities  can  be  generated  from  both  the 
baseline  FE  model  and  from  the  ABC  systems,  and  combined  in  a  global  sensitivity 
matrix  which  can  reduce  or  eliminate  conditioning  problems.  A  genetic  algorithm  (GA)- 
based  methodology,  to  make  effective  use  of  the  ABC  frequencies  for  FE  model 
updating,  is  presented  in  [14]  by  Lu  and  Tu,  where  a  binary  code,  GA,  is  proposed  for  the 
selection  of  the  desired  artificial  boundary  conditions.  The  recent  literature  related  to 
ABCs  in  general  is  completed  with  [15],  where  the  practical  aspects  of  the 
implementation  of  ABCs  in  an  actual  measurement  environment  is  investigated  and  the 
feasibility  of  using  only  the  lowest  few  ABC  frequencies  in  conjunction  with  a  variety  of 
ABC  configurations  is  examined. 

C.  REFERENCE  TO  PRESENT  WORK 

In  this  thesis,  a  sensitivity-based  FE  analysis  is  perfonned,  with  the  goal  of 
eliminating  the  discrepancies  mentioned  above.  First  of  all,  the  properties  of  several 
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Matlab  algorithms  are  examined  in  order  to  choose  the  one  that  best  matches  our 
requirements  for  linear  system  solution.  These  algorithms  are  invoked  with  the 
commands  “pinv”  and  ‘\”  and,  eventually,  involve  SVD  and  QR  decomposition  with 
column  pivoting,  respectively.  In  addition,  the  correlation  between  the  sensitivity  matrix 
rows  and  modal  strain  energy  is  depicted.  Then,  the  degree  of  the  data  incompleteness, 
which  is  related  to  the  number  of  available  modes,  is  connected  to  the  qualitative  result  of 
the  linear  system  solution.  This  leads  further  to  the  conclusion  that  there  is  a  necessary 
but  not  sufficient  condition  for  detecting  potential  damage  everywhere  in  a  structure, 
using  only  one  sensitivity  matrix  for  all  the  cases.  The  condition  is  that  one  must  provide 
as  many  modes  as  the  adjustable  design  parameters  of  the  structure,  i.e.,  building  a  square 
sensitivity  matrix.  Toward  this  end,  several  methods  of  data  enrichment  are  taken  into 
account,  like  eigenvector  sensitivities  and  clustering/grouping  parameters,  with  emphasis 
on  the  use  of  ABCs,  since  they  are  considered  the  most  effective  tool  for  all  the 
aforementioned  reasons.  Keeping  in  mind  that  even  with  a  potentially  large  number  of 
ABC  configurations  available,  a  sufficient  amount  of  additional  data  from  the  ABCs  may 
not  be  sufficient,  such  as  in  very  large  and  complex  structures;  in  these  cases,  ABC 
configurations  can  be  combined  with  grouped  parameters  into  a  global  sensitivity  matrix. 
Finally,  it  is  shown  that  the  optimal  configurations  result  in  a  diagonal  dominant  global 
sensitivity  matrix,  which  once  built,  will  meet  necessary  and  sufficient  conditions  for 
damage  detection. 
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II.  THEORY 


Beginning  from  the  governing  equation  that  describes  a  dynamic  system  under  an 
excitation,  we  have: 

[M]{x]  +  [C]{x]  +  [*]{*}  =  {F(0}  (2. 1) 

where  this  system  consists  of  mass,  damping  and  stiffness  parameters,  expressed  as  the 
symmetric  mass,  damping  and  stiffness  matrices  shown  above,  of  size  n  x  n,  where  n  is 
the  number  of  degrees  of  freedom  (DOF)  in  the  model.  The  vector  jF(t)}  represents  the 
excitation  force  of  size  n  xl. 

If  the  damping  [C]  is  neglected,  Equation  (2. 1)  becomes: 

[M]{x}  +  [^]{x}  =  {E(t)}  (2.2) 


Expressing  in  a  more  explicit  way  each  matrix  shown  in  Equation  (2.2),  results  in 
the  following: 


(2.3) 


Equation  (2.3)  describes  a  system  of  n  degrees  of  freedom,  which  approximates  a 
continuous  system  with  an  infinite  number  of  degrees  of  freedom.  An  actual  modal  test 
involves  the  measurement  of  only  a  small  percentage  of  a  structure’s  DOFs.  During  this 
procedure,  these  measured  coordinates  are  considered  the  analytical  set  or  (ASET),  and 
they  are  the  locations  where  the  measuring  instruments  (accelerometers)  are  positioned. 
The  remaining  unmeasured  DOF  are  considered  the  omitted  coordinate  set  or  (OSET), 
and  in  theory,  contain  an  infinite  number  of  degrees  of  freedom. 
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A. 


ANALYTICAL  AND  OMITTED  COORDINATE  SETS 


Going  from  the  time  domain  to  the  frequency  domain  and  getting  the  steady  state 
response  of  the  above  described  system,  we  get  for  the  whole  set  of  n  DOF: 
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where  Q  is  the  frequency  of  the  harmonic  excitation.  Making  use  of  the  analytical 
(measured)  and  the  omitted  (unmeasured)  set  of  coordinates,  that  is  the  ASET  and  OSET 
coordinates,  Equation  (2.4)  takes  the  following  form: 
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(2.5) 


Without  loss  of  any  generality,  we  can  assume  that  there  is  no  excitation  on  the  omitted 
coordinates,  that  is  {/0}  =  {0}  ,  and  get: 
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(2.6) 


The  above  matrix  form  represents  a  system  of  two  equations.  From  those  equations,  there 
is  an  obvious  solution  for  the  OSET  coordinates  in  terms  of  ASET  coordinates  through: 

[*«,  ]k  }+  [Kao  ]k  }-  Q2  [[Moa  ]{xa }+  \Moo  ]{xa }]  =  0  (2.7) 

Rearranging  and  solving  for{x0}  ,  Equation  (2.7)  becomes: 

fc }  =  [/  -  [-  K-'CK„,  +n'-K2Mm  }  (2.8) 

By  definition,  the  tenn  [y  -  Q 2Koo~lMoo  j  is  equal  to: 
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(2.9) 


I  -Cl2K  XM 

00  oo 


1 

Det[l-n2K-oloMoo\ 


Adj[l-tfK;lMoo\ 


where  the  Det[ I  -CL2Knn  xMnn ]  indicates  the  determinant  and  Adi  T/  -CL2Knn  1 

OO  OO  I  a  oo  OO  I 

indicates  the  adjoint  matrix.  From  Equation  (2.9),  it  is  obvious  that  there  are  some 
frequency  values  for  which  the  inverse  does  not  exist,  and  these  are  the  ones  that  satisfy: 

Det  [/  -  Q 2  A'J  M uo  \  =  0  (2.10) 


It  is  easily  noticeable  that  the  frequencies  that  satisfy  the  above  equation  are  the 
eigenvalues  of  the  system  defined  by  the  [ifyjand  [Ml  matrices,  hence,  they 
correspond  to  the  natural  frequencies  of  the  OSET  system. 

Realizing  that  the  eigensolution  of  a  dynamic  system  is  defined  by  the 
unrestrained  DOF  of  the  corresponding  system,  the  OSET  coordinates  are  unrestrained 
and  the  ASET  coordinates  are  considered  constrained  to  the  ground. 

B.  REDUCED  ORDER  MODEL 


The  fact  that  a  finite  number  of  transducers  are  used  for  a  modal  test,  in  contrast 
to  the  infinite  number  of  DOF  of  the  real  system,  defines  a  reduced  order  model.  This 
model’s  impedance  is  nonlinearly  dependent  on  the  impedance  of  the  full  order  model 
[13].  Considering  this  exact  model,  in  tenns  of  frequency  response  function  (FRF)  we 
have: 
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(2.11) 


The  experimental  set  of  data  measured  in  a  real  test  is  carried  in  the  partition 
\Haa]  of  the  above  full  FRF.  In  other  words,  [//M]  represents  a  structural  dynamic  model, 

reduced  via  exact  dynamic  reduction  [11].  Using  the  algebraic  relation  of  FRF  to  its 
impedance,  we  have: 
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(2.12) 
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which,  after  decomposing  the  set  of  equations  and  solving  for  \Haa  ] ,  leads  to: 


[vj=[z 


(2.13) 


Writing  Equation  (2.6)  in  terms  of  impedance,  we  get: 


(2.14) 


Rearranging  (2.14),  the  exact  dynamic  relationship  between  the  OSET  coordinates  and 
the  ASET  coordinates  gives: 


(2.15) 


Combining  Equation  (2.14)  and  Equation  (2.15),  and  solving  the  first  equation  of  the 
matrix  formula,  the  dynamic  reduced  model  results: 

{/  }  =  [z  -Z  Z_1Z  l{x  I  (2.16) 

t  J  a  )  I  aa  ao  oo  oa  J  (  a )  v  z 

Paralleling  Equations  (2.13)  and  (2.16),  it  is  found  that  the  FRF  extracted  from  an 
actual  modal  test,  i.e.,  the  measured  partition  of  an  infinite-order  FRF  matrix,  is 
equivalent  to  the  inverse  of  a  dynamically  reduced  impedance  matrix.  In  addition,  the 
term  Zoo  1  seems  to  be  critical  for  the  definition  of \Haa  1 ,  and  since  every  element  in  the 

Zoo  '  is  singular  at  the  natural  frequencies  of  the  OSET,  Equation  (2.13)  shows  that  the 
elements  of  Haal  will  be  singular  at  the  OSET  natural  frequencies  [13]. 
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From  what  is  proved  above,  it  is  obvious  that  the  FE  model  can  be  manipulated  in 
the  same  manner.  This  means  that  the  FE  model  is  used  to  calculate  a  full  n  x  n  FRF,  and 
appropriate  boundary  conditions  can  be  applied  corresponding  to  the  coordinates 
measured  in  a  test. 

C.  FREQUENCY  RESPONSE  FUNCTION 

The  Frequency  Response  Function  FRF  is  actually  a  ratio  defined  by  the  response 
of  a  system  due  to  a  harmonic  excitation.  It  contains  all  the  amplitude  and  phase 
information  for  the  response  of  a  specific  DOF  excited  by  an  input  force  at  a  specific,  and 
perhaps  different,  DOF.  The  FRF  is  constructed  by  measuring  both  the  response  and  the 
excitation  simultaneously.  This  data  is  transformed  from  the  time  domain  to  the 
frequency  domain  using  FFT  algorithm.  Doing  this  for  all  the  DOF  that  are  set  to  be 
measured  for  certain  external  excitations,  the  FRF  matrix  [//(Q)]  can  be  built.  As 

stated  before,  the  inverse  of  this  matrix  is  known  as  the  impedance  matrix  [Z  (Q)]  : 

[Z(fi)]  =  (2.17) 


where: 


[Z(Q)]  =  [K-DrM  +  jO.C] 


(2.18) 


In  tenns  of  an  n  DOF  system  and  introducing  impedance,  Equation  (2.4)  can  be  written 
as: 
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(2.19) 


Letting  jx}  =  [O]  (c/j  ,  where  O  is  the  mass  normalized  mode  shapes  of  the  system  and 
<7  represents  the  generalized  coordinate  set,  substituting  in  (2.19): 

FM?1={/} 
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(2.20) 


Then  premultiplying  by  of  and  combining  Equations  (2.18)  and  (2.20): 


M  M$] -Q2  [of  [M][o]  +  jD[of  [C][0]  {q}  =  [Of  {F}  (2.21) 


Since  O  is  mass  normalized  [of  [M][0]  =  1 ,  and  assuming  proportional  damping 
[C\  =  a*[K]  +  J3*[M]: 


\co?-D2  +  2 =  [of  {F}  (2.22) 


In  Equation  2.22,  coi  is  the  i,h  natural  frequency,  D  is  the  forcing  frequency  and 
C,  is  the  systems  damping  ratio.  The  term  'far  -D2  +  I  jCDojf  is  known  as  the  modal 

impedance  matrix  and  is  diagonal.  Inverting  the  modal  impedance  matrix  and  moving 
back  to  physical  coordinates  by  premultiplying  by[0] ,  the  frequency  response  is: 


=M 


a>?  -D2  +  2j<^Dcof 


[®f{F} 


(2.23) 


From  Equation  (2.23),  the  FRF  matrix  [//(Q)]  can  be  defined  as: 


[//(Q)]=[0] 


co?  -Q2  +  2  j<^Dcoi 


N7 


(2.24) 


The  FRF  matrix  can  be  rewritten  in  summation  form: 

\ok\ok 


modes 


[»(«)]=  £ 


k=l  ^  ^  J  $^^k 


(2.25) 


and  descriptively  for  any  element  in  terms  of  Htj  (Q)J  : 
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(2.26) 


kN-z  t tSL 

k= i  ~  O  +  2j£Q.0Jk 

It  is  obvious  from  Equation  (2.26)  that  the  FRF  contains  all  the  modal  data  for  the 
system. 


1.  Driving  Point  and  Transfer  Point  Functions 

The  exact  FRF  matrix,  as  mentioned,  contains  all  the  modal  data  of  the  system 
and  thus  the  resonance  and  anti-resonance  frequencies,  which  are  of  particular  interest  in 
modal  analysis.  It  should  be  noted  here  that,  depending  on  the  relationship,  in  terms  of 
locality,  of  the  ith  DOF  where  the  response  is  measured  with  the  j"‘  DOF  where  the  force 
is  applied,  the  FRF  is  referred  to  as  the  Driving  Point  Function  for  the  case  that  these  two 
DOF  are  identical  and  Transfer  Function  for  all  the  other  cases. 

What  follows  is  an  example  of  FRF  derivation,  covering  both  Driving  Point  and 
Transfer  Function,  for  a  simple  two  DOF  system  with  no  damping  [16]. 


Figure  1.  2-DOF  Example 

The  properties  of  this  system  are  summarized  in  the  stiffness  and  mass  matrices 

above: 
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Let,  for  this  example,  Mx  =  M,  =1.0  and  Kx  =  K3  =0.4  and  K2  =0.8.  This  yields  the 
following  results: 
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Mode  shapes:  |<E>'| 


J— 0.707  ll 
[-0.707  lj 


and 


[-0.7071] 
[0.7071  J 


Natural  frequencies:  [co(rad  /  sec)} 


[0. 6325} 
[1.0954  J 


Using  these  values  in  Equation  (2.26),  //^.(Q)] 


h  col  -  D.1  +  2  jCD.0)k 


results  in: 


M{< ,  MKT 

COI-  O2  co;-n2 


{-0.7071}  {-0.707l}r  {-0.7071}  {-0.707l}r 

0.4  -Q2  1.2  -Q2 


(2.27) 


for  the  Driving  Point  function  and  for  the  Transfer  function  in: 


[//12(Q)] 


{o;}{o[}r  {o{}{0 
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col -n2  col-n2 


{-0.7071}  {-0.707l}r  {-0.7071}  {0.707l}r 

0.4 -Q2  1.2  -Q2 


(2.28) 


What  is  evitable  from  Equation  (2.27)  for  the  driving  point  is  that  for  frequencies 
below  and  above  the  natural  frequencies  of  the  system,  the  two  terms  will  have  the  same 
sign,  while  between  them,  they  will  have  the  opposite  sign.  This  means  that  in  the  first 
case,  the  two  terms  are  additive,  where  in  the  second  they  are  subtractive.  Furthennore, 
this  effect  creates  resonances  at  forcing  frequencies  equal  to  the  natural  frequencies  and 
for  the  particular  values  between  them,  the  two  terms  of  the  equation  may  result  in  zero 
value,  thus  points  of  anti-resonance  are  created.  Therefore,  anti-resonance  follows  every 
resonance,  without  exception,  as  seen  in  the  top  portion  of  Figure  2. 

On  the  other  hand,  there  is  no  reason  to  expect  something  similar  for  the  Transfer 
Function,  because  obviously  the  two  terms  of  Equation  (2.28)  will  sometimes  be  additive 
and  sometimes  be  subtractive,  as  the  product  of  the  two  eigenvectors — called  modal 
constant — will  be  sometimes  positive  and  sometimes  negative.  Thus,  we  expect  Transfer 
FRF  measurements  to  show  a  mixture  of  anti-resonances  and  minima.  It  can  be  shown 
that,  in  general,  the  further  apart  the  two  points  in  question  are,  the  more  likely  the  two 
eigenvector  elements  are  to  alternate  in  sign  as  one  progresses  through  the  modes  [17]. 
The  above  conclusions  are  shown  in  the  lower  portion  of  Figure  2. 
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Driving  Point  FRF,  H11 


Transfer  FRF,  HI 2 


Figure  2.  2-DOF  Frequency  Response  Function. 

The  principles  demonstrated  in  the  above  example  can  be  applied  to  any  number 
of  DOF.  This  concept  of  the  existence  of  anti-resonance  between  any  two  modes  of  a 
driving  point  FRF  is  of  great  significance  in  the  use  of  Artificial  Boundary  Conditions. 

D.  ARTIFICIAL  BOUNDARY  CONDITIONS 

As  discussed  previously,  the  generation  of  additional  data  for  a  spatially 

incomplete  FE  model  analysis  can  be  pursued  via  Artificial  Boundary  Conditions,  or 

ABC.  This  way,  a  number  of  ABC  can  be  applied  to  an  FE  model,  thus  constructing 

several  models  differing  only  in  boundary  conditions.  Hence,  a  new  model  configuration 

is  provided,  yet  only  one  set  of  measured  data  is  required.  Applying  the  driving  point 

function  concept  shows  that  the  ABC  frequencies  of  a  certain  ASET  coordinate  are 

defined  by  the  corresponding  driving  point  anti-resonances.  In  other  words,  these  ABCs 

are  the  boundary  conditions  that  define  the  OSET,  from  which  it  is  shown  that  additional 

frequency  information  about  the  model  is  provided.  The  following  examples  describe 
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how  anti-resonance  frequencies  are  identical  to  the  natural  frequencies  corresponding  to 
the  system  with  the  associate  driving  point  constrained  to  the  ground  [13]. 

1.  ABC  Frequencies  for  a  2-DOF  System 

Consider  the  system  shown  in  Figure  3,  excluding  damping: 


Figure  3.  2-DOF  System 


From  Equation  (2.26),  the  driving  point  FRF  is 


modes 


k= 1 


($-)2 

co2k -Or 


(2.29) 


where  dy  is  the  mass  normalized  mode  shape  element,  co*  is  the  &th  natural  frequency, 
and  Q  is  the  forcing  frequency. 

Solving  for  the  frequency  of  the  anti-resonance  of  Hu(  Q)  yields: 
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2 

anti-res 


R\\Col  +  R^cof 
*n+*n 


(2.30) 


where  the  modal  residue  is  given  by  Ry  =  d>*d>^ . 

Let  Ml  =  M2  =  1.0  and  K1  =  K2  =1.0.  Then  the  stiffness  and  mass  matrices  of 
the  system  are: 
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Solving  the  eigenproblem  related  to  them  yields  to: 
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Natural  frequencies:  {a>(rad  /  sec)} 


JO. 6181 
[1.618  J 


and  the  single  anti-resonance  at  Clantires  =  a/2  rad/sec,  as  shown  in  Figure  4. 


Driving  Point  FRF,  H11 
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Figure  4.  Driving  Point  FRF,  HI  1 
Now,  consider  the  same  system  pinned  at  DOF  1,  as  shown  in  Figure  5: 


15 


Figure  5.  DOF  #1  Constrained  to  Ground 


This  system’s  natural  frequency  is  co  =  yj 2  rad/sec  and  identically  equal  to  the 
anti-resonance  of  the  unconstrained  system  of  Figure  3.  The  same  result  for  the 
constrained  system  can  be  found  by  making  use  of  Equation  (2.13),  that  is,  by  calculating 
the  inverse  of  the  ASET  FRF  partition  that  is  constrained//”'  (Q) ,  which  in  this  case  is 

the  inverse  of//n(Q).  The  correlation  between  the  driving  point  anti-resonance  and  the 
natural  frequency  of  the  system  with  the  driving  point  constrained  to  ground  is  shown  in 
Figure  6  [13]. 


Plot  A.  Driving  Point  FRF.  H11 


Plot  B.  Z  =inv(H1 1 )  of  the  ABC  System 


Figure  6.  Plot  A.  Driving  Point  HI  1(D)  of  System  1,  Plot  B.  [Haa(D)]-l  of  System  2 
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2. 


ABC  Frequencies  for  a  Free-Free  Beam 


The  model  depicted  in  Figure  7  consists  of  ten  beam  elements  and  each  element 
contains  two  nodes  with  two  DOF  each,  the  odd-numbered  DOF  being  translational,  and 
the  even  numbered  DOF  being  rotational.  Response  transducers  are  set  up  at  the 
translational  DOF:  1,  5,  9,  13,  17,  and  21  as  follows: 


_ 
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71 
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8 
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10 

Figure  7.  Transducers  Located  at  DOF  1,  5,  9,  13,  17,  and  21 

The  location  of  the  response  accelerators  represent  the  ASET,  defined  as  [1,  5,  9, 
13,  17,  21].  If  it  is  assumed  that  the  excitation  is  applied  at  each  of  the  ASET  DOF,  then  a 
square  (6  x  6)  FRF  matrix  will  result  from  the  testing.  Then  the  impedance  matrix 
H~la{Q.)  is  calculated  in  the  bandwidth  0-800  Hz.  Figure  8  shows  the  driving  point  FRF 
of  the  system. 


Frequency  (Hz) 


Figure  8.  Driving  Point  FRF  ASET  DOF:  1  5  9  13  17  21.  (From  [13]) 
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Now,  consider  the  same  system  with  all  the  ASET  DOF  constrained  to  the 
ground,  as  shown  in  Figure  9. 


Figure  9.  ABC  Configuration  ASET  [1  59  13  1721]  (Measured  Coordinates  Restrained 

to  Ground) 


By  calculating  the  impedance  matrix  7/J(G)  and  plotting  the  //, ,  (Q)  of  the  matrix  in  the 
same  bandwidth  of  0-800  Hz,  Figure  10  is  generated: 


Frequency  (Hz) 

Figure  10.  Impedance  FRF  for  ASET  DOF:  1  5  9  13  17  21  (From  [13]) 

By  deriving  the  ABC  frequencies  of  the  system  in  Figure  9,  it  is  obvious  that  they 
correspond  exactly  to  the  resonances  of  the  impedance  FRF  for  ASET  coordinates  the 
DOF,  which  are  constrained  to  the  ground. 
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E. 


SENSITIVITY  BASED  FEM  UPDATING 


Sensitivity  analysis  in  structural  dynamics  is  a  tool  used  for  predicting  the  effect 
of  small  changes  in  design  parameters  on  dynamic  response  characteristics.  It  involves 
finding  how  a  system  alters  its  dynamic  response  when  some  parameters  that  define  its 
properties  are  subject  to  a  change.  From  a  mathematical  point  of  view,  it  implies 
differentiation,  which  in  the  present  analysis,  is  carried  out  through  finite  difference 
calculations.  The  parameters  of  the  system  to  be  changed,  stated  as  design  variables,  are 
in  this  case,  only  stiffness  parameters,  whereas  other  parameters,  like  mass  or  density  can 
be  taken  into  consideration.  The  output  of  the  system  to  be  examined  using  sensitivities, 
due  to  an  input  parameter  change,  can  be  one  or  both  of  the  following:  natural 
frequencies,  eigenvalue  sensitivity  and  mode  shapes,  eigenvector  sensitivity.  Another 
aspect  of  modal  data  analysis  via  sensitivity  is  the  Frequency  Response  Function 
Sensitivity  Analysis,  which  is  not  covered  in  this  thesis. 

Recalling  the  ABC  concept,  it  is  shown  that  ABC  configurations  can  be  used 
along  with  the  baseline  system  into  a  combined  sensitivity  matrix,  for  model  updating 
[13].  The  key  part  of  this  is  that  ABC  frequencies  can  easily  be  found,  as  stated 
previously,  and  since  the  design  variables  are  common  for  both  the  ABC  and  the  baseline 
system,  they  can  be  used  to  reliably  enlarge  the  modal  data  for  the  analysis. 

The  equation  that  describes  what  is  mentioned  above  about  sensitivity  based 
updating  is: 

{A«2}=[r]{ADC}  (2.31) 


where  |a<u2|  is  a  vector  containing  the  difference  between  the  experimental  and  the 
analytical  (FE  baseline)  natural  frequencies  jo2  -®2}  •  The  (ADFj  term  is  the  vector 
containing  the  design  variables,  and  [T ]  is  the  sensitivity  matrix. 

Each  column  of  the  sensitivity  matrix  represents  an  element  of  the  model, 
whereas  each  row  represents  a  mode  of  the  model. 
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Each  ABC  system  defines  additional  rows  to  Equation  (2.31).  It  is  obvious  from 
this  point  that  even  with  low  modes,  in  terms  of  frequency  spectrum,  that  is  with  the  first 
few  modes  of  baseline  and  ABC  systems,  a  more  determined  sensitivity  matrix  can  be 
built.  Representing  the  baseline  system  quantities  with  the  superscript  “0”  and  the  ABC 
systems  with  the  superscripts  from  1  to  “k,”  where  k  represents  the  number  of  ABC  used 
in  the  system,  the  following  equation  shows  how  this  is  compromised  [13]: 


A 

( T’(O)  ^ 

1  BASE 

r 

•33  ... 

A<4ci 

>  = 

T’i  1) 

1  ABCl 

< 

A<^ 

ji(k) 

K1  ABCk ) 

m 

(2.32) 


1.  Eigenvalue  Sensitivity  Matrix  Derivation. 

The  starting  point  of  the  derivation  is  the  eigenvalue  problem  of  an  n  DOF 

system: 


(2.33) 


Differentiating  with  respect  to  the  design  variable  {DV}: 

A Ai 


A  K  AM _ 

A DV  'A DV  A DV 


(2.34) 


Expanding  Equation  (2.34)  and  premultiplying  by  {0}r  leads  to: 


“  A K  " 

“  AM  " 

“  A Ai  " 

{<3>,}  -  A,  {<!>,} 

{O’,}- 

A DV  _ 

A DV 

t  ) 

A DV 

=  (2.35) 


By  invoking  the  matrix  identity  property  {a}r  [£>]  {c}  ={c}r[b]{a},  the  last  term  on  the 
left  hand  side  can  be  replaced  by: 


AO,  Y 
ADV} 


[K-hM]{  O,}  =  0 


(2.36) 
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Since  [K  -/IjWjjOj  =  {0}  and  recalling  the  mass  normalization  of  the  mode  shapes, 
where  [®;]7  [M][®,]  =  1  ,  the  overall  equation  is  reduced  to: 


®, 


A K 
A DV 


{®,}-/k{®i}7 


AM 
A  DV 


{«>'}- 


AAa 
A DV 


=  M 


(2.37) 


Bringing  the  mass  and  stiffness  portion  to  the  right  hand  side 


“  A  Ai  " 

(  ^  \T 

“  A K 

.  AM 

={®<} 

A DV  _ 

A DV 

A DV 

(2.38) 


As  stated  previously,  the  present  analysis  only  takes  into  account  changes  in  the 

stiffness  parameters.  Therefore,  the  second  term  h  in  the  square  brackets  of  the 

A DV 

right  hand  side  can  be  neglected  as  AM  =  0 ,  yielding  to: 

{®/}  where  [AA]  =  [Ah- A«]  (2.39) 


A/k 
A DV 


=  1®, 


A  K 
A DV 


and  finally: 


A K 
A DV 


{®4 


(2.40) 


It  has  been  shown  [18]  that  when  ABC  are  applied,  the  model  tends  to  concentrate 
high  stiffness  sensitivity  values  around  the  particular  DOF  where  the  ABC  is  located. 
This  concept,  along  with  the  fact  that  the  baseline  model  is  not  likely  to  predict  changes 
in  natural  frequencies  in  regions  with  low  sensitivity  values,  is  of  great  significance  in  the 
analysis  of  the  following  chapters. 

2.  Eigenvector  Sensitivity  Matrix  Derivation 

The  calculation  of  the  eigenvector  sensitivity  is  a  less  straightforward  problem 
than  that  of  the  eigenvalue  sensitivity.  The  foundation  for  the  calculation  was  laid 
primarily  by  Fox  and  Kapoor  [19]  and  a  simplified  method,  which  is  the  basis  of  the 
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NASTRAN  computation  scheme,  is  provided  later  by  Nelson  [20].  A  survey  of  several 
contributions  to  the  literature  regarding  the  eigenvector  sensitivity  calculation  is  available 
in  [21]. 

The  derivation  that  follows  is  based  on  Blelloch’s  report  [22].  After  having 
calculated  the  eigenvalue  sensitivity  from  Equation  (2.39),  returning  to  Equation  (2.34) 
and  solving  for  the  eigenvector  derivative,  again  for  no  mass  variation  yields  to: 


A K 


A  Xi 


A DV  A DV 


M 


(2.41) 


Assume  at  this  point  that  all  the  eigenvectors  are  available  for  the  system  under 
consideration.  In  this  case,  the  eigenvector  derivative  can  be  expressed  as  a  linear 
combination  of  all  the  eigenvectors: 


[AAl 

A DV 


=  [<■>]{< 


(2.42) 


where  [<f>]  is  the  matrix  containing  all  the  eigenvectors  and  jc}  is  the  vector  of  the 
coefficients  to  be  determined.  Substituting  to  Equation  (2.41)  and  premultiplying  by 
M  generates: 


A  K 
A DV 


A  DV 


Ml 


(2.43) 


The  expression  that  premultiplies  jc}  is  of  the  form: 

K[a-  /L/M][<h]  =  [A  -  /IT] 


(2.44) 


where  A  is  the  diagonal  matrix  of  eigenvalues.  The  above  implies  that,  for  distinct 
eigenvalues,  every  element  of  the  coefficient  vector  jc}  except  for  the  ith  element  ( c  )  is 
uniquely  detennined  as  follows: 
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c 
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A K  A  h  , , 

- M 

ADV  ADV 

\  ~  ^ i 


(2.45) 


At  this  point,  it  should  be  noted  that  there  is  no  reason  to  believe  that  the 
eigenvector  derivative  will  predict  the  change  in  the  eigenvector  in  models  with  closely 
spaced  eigenvalues,  since  from  what  can  be  seen  from  Equation  (2.45)  numerical  ill- 
conditioning  will  occur  in  these  cases. 

The  eigenvector  derivative  can  now  be  written  as: 

■  ^7} = to) + Ml  =  to! to )  <2-4«) 


In  this  expression,  the  vector  |jto}  is  detennined  uniquely.  The  undetennined  coefficient 
c;  will  be  found  on  the  basis  of  differentiating  the  mass  normalization  fonnula  for  the 
eigenvectors: 

^(tofMtoH) 


to  have: 


8M 

dDV 


M}  =  0 


(2.47) 


Again,  for  no  mass  variation  and  substituting  Equation  (2.46)  into  Equation  (2.47): 

(2.48) 

The  key  to  the  whole  method  that  Nelson  [20]  indicated  is  that  there  is  no  need  to 
calculate  the  unique  {k} .  In  fact,  any  vector  that  satisfies  Equation  (2.41)  can  be  used  and 
an  available  check,  as  stated  in  [22],  is  that  the  eigenvector  derivative  is  orthogonal  to  the 
eigenvector.  Following,  in  trying  to  find  an  arbitrary  non-  trivial  solution  for  { 17} ,  the  any 
kth  row  and  column  are  partitioned  from  the  eigenproblem  matrix  as  follows: 


23 


(2.49) 


(K-X,M)U  ( K-X,M)U  (K~AM\, 

(K-A,M)n  (K-^U)n  < 

_(K-Z,M)n  (K-A,M\k  (K-A,M)3J_ 

Choosing  k  such  that  xk  is  non-zero  and  solving  for  x,  and  x3  in  terms  of  xk : 
~(K-A,M)u(K-A,M)n 

(k-a,mUk-a,m)k 

Since  the  kth  row  (equal  to  the  kth  column  due  to  symmetric  matrices)  is  linearly 
dependent  on  the  ( n-1 )  remaining  rows,  the  kth  row  can  be  eliminated  and  a  system  of  (n- 
1 )  equations  of  rank  {n-1)  is  guaranteed  as  follows: 

~(K-A,M)u  (k-a,m\, 

(K-A,M)3l  (K-A,M)n 
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■  =  -Xt 


\(K-W\ 


(2.51) 
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(2.50) 
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Evoking  that  [FJ }  is  arbitrary;  its  kth  element  can  be  set  to  zero  and  using  Equation  (2.5 1) 
to  solve  for  ( E  j  results: 
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The  unique  solution  for  the  eigenvector  derivative  will  then  be: 


J  i 
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23  J 


(2.52) 
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(2.53) 
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where  {Fj}  and  { V3 j  are  calculated  from  Equation  (2.52)  and  ci  from  Equation  (2.48). 


An  algorithm  for  calculating  the  eigenvector  derivative  (sensitivity),  assuming  eigenvalue 
sensitivity  has  already  been  calculated,  is  given  below  [22]: 

•  Loop  over  all  eigenvalues 

•  Calculate  the  left  hand  side  matrix  [Z).  ]  =  [X  -  A M] 

•  Select  k  based  on  the  largest  element  of 

•  Partition  out  kth  row  and  column  from  [  Z).  ] 


\D] 

nv  l  'J 


•  Decompose  the  matrix 

•  Loop  overall  design  variables 


A K  A  Ai 


M 


Calculate  the  RHS  vector  {  F  \  = 

1  J  1 ADV  A DV 

Forward  Backwards  Substitution  to  calculate  {if} 

Merge  {if}  to  an  //-vector  by  adding  zero  kth  element 

Calculate  ct  =  -{$}T  [M ]{F]} 

Eigenvector  derivative  is  given  by  { E }  +  c;  { i/;  j 

•  End  of  design  variable  loop 

End  of  eigenvalue  loop 


F.  CORRELATION  OF  STRAIN  ENERGY  TO  SENSITIVITY  MATRIX 


At  this  point,  it  would  be  useful  to  make  a  brief  analysis  of  the  strain  energy  of  a 
dynamic  system  and  to  examine  how  it  is  correlated  to  the  sensitivity  matrix.  A  relevant 
analysis  has  been  perfonned  in  [18],  where  the  shape  functions  of  an  Euler-  Bernoulli 
beam  were  used  for  the  derivation  of  the  strain  energy  of  the  system.  In  this  work,  after 
having  set  the  FE  model  for  the  structure,  which  requires  the  calculation  of  the  system’s 
mode  shapes  and  of  each  element  stiffness  matrix,  a  simpler  approach  was  made  using 
the  trivial  equation: 
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U=-kx 2 
2 


(2.54) 


which  for  each  element  of  the  model  using  modal  data  is  translated  into: 

Ue  =  ^-[o£']T  [if£'][(De]  (2.55) 

where  stands  for  the  element  stiffness  matrix  and  [d>e]  stands  for  the  partition  of 

the  mode  shape  matrix  that  contains  the  eigenvectors  related  to  the  DOF  associated  with 
the  particular  element.  It  is  of  vital  importance  to  be  noted  here,  that  this  mode  shape 
matrix  is  not  the  original  eigenvector  matrix  of  the  baseline  system,  in  the  way  that  it  has 
been  modified  in  order  to  take  into  account  the  physical  boundary  conditions  of  the 
system  which  may  occur  in  some  of  the  elements’  boundaries.  For  example,  consider  a 
cantilever  beam  with  its  left  end  clamped  is  examined  and  the  [d^]  corresponding  to  the 

first  element  of  the  FE  model.  The  entries  of  the  matrix  associated  with  the  right  two 
DOF  of  the  element  will  be  covered  by  the  entries  of  the  baseline  system  eigenvector 
matrix  associated  with  the  first  two  (unconstrained  as  trivially  defined)  DOF  of  the 
cantilever  beam,  while  the  entries  associated  with  the  left  two  DOF  of  the  element  must 
be  fulfilled  with  zeros  corresponding  to  the  constraint  imposed  to  both  the  translational 
and  the  rotational  DOF. 

That  said,  and  having  derived  the  stiffness  sensitivity  matrix  for  the  same  system, 
a  very  remarkable  conclusion  is  reached.  As  shown  in  the  following  Figures  11  and  12, 
where  all  the  modes  of  a  10  element  cantilever  beam  are  depicted,  there  is  an  exact  and 
direct  correlation  between  the  stiffness  sensitivity  matrix  and  the  flexural  properties  of 
the  structure.  An  appropriate  scaling  with  nonnalization  to  the  maximum  value  of  its 
mode  has  been  made  to  both  of  them  in  order  to  make  the  correlation  more  legible.  In 
these  figures,  the  x-axis  has  the  number  of  elements;  the  y-axis  has  the  normalized  value 
of  both  the  sensitivity  matrix  and  the  strain  energy,  which  are  referred  to  each  mode  of 
the  system,  respectively. 
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Correlation  of  Sensitivity  Matrix  to  Strain  Energy  for  Modes  1  Through  10 
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Figure  12.  Correlation  of  Sensitivity  Matrix  to  Strain  Energy  for  Modes  1 1  Through  20 
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III.  LINEAR  SYSTEM  SOLUTIONS  IN  SENSITIVITY  BASED 

FEM  UPDATING 


As  demonstrated  in  the  previous  chapters,  the  procedure  of  sensitivity-based  FEM 
updating  involved  the  solution  of  a  linear  system,  which  was  almost  always 
underdetennined.  This  linear  system  is  governed  by  the  Equation  (2.31): 

{Ary2}=[r]{ADE} 

In  FEM  updating,  the  vector  {ADF}  contains  corrections  to  the  model 
parameters,  while  in  damage  detection,  the  vector  { AD  V}  contains  flags  to  both  the 
location  and  severity  of  structural  damage,  relative  to  the  undamaged  FEM.  In  both 
cases,  Equation  (2.31)  has  to  be  solved  for  the  design  variables  { AD V] .  This  process 
constitutes  an  inverse  type  of  mathematical  problem,  which  requires  the  inversion  of  the 
sensitivity  matrix \T ] .  Of  course,  this  is  feasible  only  if  the  matrix  is  square  and  full  rank, 

which  is  rarely  the  case.  For  all  cases  where  the  sensitivity  matrix  is  underdetennined, 
there  are  several  algebraic  methods  of  linear  systems  for  solving  the  problem,  which  are 
covered  in  the  literature  [2], 

A.  USING  BACKSLASH  “\”  AND  “PINV”  IN  MATLAB 

In  this  work,  Matlab  is  used  to  carry  out  this  procedure.  In  dealing  with  these 
kinds  of  systems,  there  are  two  available  commands/operators  that  can  provide  the 
required  pseudo-inverse  of  the  sensitivity  matrix;  these  are  “\”  or  “mldivide,”  and  “pinv.” 
Respectively,  they  use  orthogonal-triangular  (QR)  decomposition,  and  singular  value 
decomposition  (SVD).  Eventually,  they  both  give  a  least-squares  (square  and  full  rank 
overdetermined  cases)  or  a  least-norm  (underdetennined  and  rank  deficient  cases) 
solution  to  Equation  (2.3 1),  which  is  pursued  by  minimizing  the  norm  of  the  residual,  i.e., 

min|[r]{ADF}-|A(U2|  ,  except  that  they  do  it  in  different  ways.  From  what  has 
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already  been  discussed,  it  is  obvious  that  for  rank  deficient  problems,  the  solution 
obtained  from  each  one  is  not  exact.  Some  brief  descriptions  of  how  each  of  them  is  used 
[23]  are  given  below. 

1.  Using  “\”  or  “mldivide 

The  mldivide(A,B)  and  the  equivalent  A\B,  perform  matrix  left  division  (back 
slash).  A  and  B  must  be  matrices  that  have  the  same  number  of  rows,  unless  A  is  a  scalar, 
in  which  case  A\B  performs  element-wise  division.  If  A  is  a  square  matrix,  A\B  is  the 
same  as  inv(A)*B,  where  inv(A)  stands  for  the  inverse  of  A,  except  it  is  computed  in  a 
different  way.  If  A  is  an  n-by-n  matrix  and  B  is  a  column  vector  with  n  elements,  or  a 
matrix  with  several  such  columns,  then  X  =  A\B  is  the  solution  to  the  Equation  AX  =  B 
computed  by  Gaussian  elimination  with  partial  pivoting.  A  warning  message  is  displayed 
if  A  is  badly  scaled  or  nearly  singular.  If  A  is  an  m-by-n  matrix  and  B  is  a  column  vector 
with  m  components,  or  a  matrix  with  several  such  columns,  then  X  =  A\B  is  the  solution 
in  the  least  squares  sense  to  the  overdetermined  system  of  Equations  AX  =  B,  only  if  the 
rank  k  is  equal  to  n  (number  of  columns).  In  other  words,  X  minimizes  norm(A*X  -  B), 
the  length  of  the  vector  AX  -  B.  The  rank  k  of  A  is  determined  from  the  QR 
decomposition  with  column  pivoting.  If  k  <  n  or  the  system  is  underdetermined,  this  is 
not  the  same  solution  as  x  =  pinv(A)*B,  which  returns  a  least-norm  solution,  but  in  this 
case  X=A\B  returns  a  solution  that  is  called  “basic”  solution  and  has  at  most  k  non-zero 
elements  per  column. 

2.  Using  “pinv” 

The  syntax  B  =  pinv(A)  returns  the  Moore-Penrose  pseudoinverse  of  A.  The 
Moore-Penrose  pseudoinverse  is  a  matrix  B  of  the  same  dimensions  as  A’  (transpose  of 
A)  satisfying  the  four  following  conditions:  A*B*A  =  A,  B*A*B  =  B,  A*B  and  B*A  are 
Hermitian.  The  computation  is  based  on  singular  value  decomposition  of  A.  If  A  is 
square  and  not  singular,  then  pinv(A)  is  an  expensive  way  to  compute  inv(A).  If  A  is  not 
square,  or  is  square  and  singular,  then  inv(A)  does  not  exist.  In  these  cases,  pinv(A)  has 
some,  but  not  all  of  the  properties  of  inv(A).  If  A  has  more  rows  than  columns  and  is  not 
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of  full  rank,  then  the  overdetermined  least  squares  problem  of  minimizing  norm(A*x-b) 
does  not  have  a  unique  solution.  Two  of  the  infinitely  many  solutions  are  x  =  pinv(A)*b 
and  y  =  A\b.  These  two  are  distinguished  by  the  facts  that  norm(x)  is  smaller  than  the 
norm  of  any  other  solution  and  that  y  has  the  fewest  possible  non-zero  components  as 
described  for  the  use  of 

3.  Example  Demonstrating  the  Difference  Between  the  Two  Commands 

The  following  example  is  totally  representative  of  the  properties  of  the  two 
Matlab  operators  for  this  kind  of  linear  systems  solving  [23]. 

Assume  the  matrix  generated  by  A  =  magic(8);  A  =  A(:,l:6).  This  is  an  8-by-6 
matrix  that  happens  to  have  rank(A)  =  3: 


64 

■  2 

3 

61 

60 

6 

'260 

9 

55 

54 

12 

13 

51 

260 

17 

47 

46 

20 

21 

43 

260 

40 

26 

27 

37 

36 

30 

,  and  the  RHS  is  b  =  260*ones(8,l),  that  is :b  =  < 

260 

32 

34 

35 

29 

28 

38 

260 

41 

23 

22 

44 

45 

19 

260 

49 

15 

14 

52 

53 

11 

260 

8 

58 

59 

5 

4 

62 

260 

The  scale  factor  260  is  the  8-by-8  magic  sum.  With  all  eight  columns,  one 
solution  to  A*x  =  b  would  be  a  vector  of  all  Is.  With  only  six  columns,  the  equations  are 
still  consistent,  so  a  solution  exists,  but  it  is  not  all  Is.  Since  the  matrix  is  rank  deficient, 
there  are  an  infinite  number  of  solutions.  Two  of  them  are  x  =  pinv(A)*b,  which  is: 


x  =  < 


1.1538' 

4 

1.4615 

5 

1.3846 

0 

1.3846 

> ,  and  y  =  A\b  which  produces  the  result:  y  =  < 

0 

1.4615 

0 

1.1538 

-1 

with  the  following 


warning:  Rank  deficient,  rank=  3,  tol=  1.8829e-013. 
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Both  of  these  are  exact  solutions  in  the  sense  that  norm(A*x-b)  and  norm(A*y-b) 
are  on  the  order  of  round-off  error.  The  solution  x  is  special  because  norm(x)  =  3.2817  is 
smaller  than  the  norm  of  any  other  solution,  including  nonn(y)  =  6.4807.  On  the  other 
hand,  the  solution  y  is  special  because  it  has  only  three  non-zero  components. 

From  the  statement  above,  it  is  clear  that  in  terms  of  damage  detection  and  error 
localization,  the  most  preferable  operator  for  solving  for  the  linear  system  corresponding 
to  the  structure’s  model  is  the  “\”,  as  it  returns  the  minimum  number  of  non-zero 
components.  Using  only  one  mode  that  is  one  sensitivity  matrix  row,  it  will  return  only 
one  non-zero  element,  which  for  the  right  choice  of  mode,  will  be  the  element  with  the 
true  error  (as  will  be  shown  later,  the  right  mode  is  the  one  whose  highest  value  happens 
to  be  at  the  element  where  the  true  error  occurs);  whereas  “pinv”  will  return  a  vector 
scaled  by  the  mode  used.  Using  more  than  one  mode  will  provide  the  user  with  a  more 
transparent  result  in  terms  of  localization  until  a  square  matrix  is  reached  where  there  is 
no  difference  between  the  two  commands,  provided  the  matrix  is  not  rank  deficient.  This 
may  prove  to  be  of  vital  importance  during  the  pursuing  of  optimal  configurations, 
especially  because  it  makes  clear  that  modes  of  systems/configurations  with  locally 
concentrated  sensitivity  values  will  reveal  errors  at  the  corresponding  elements  where 
these  values  are  laid. 

B.  SPATIALLY  INCOMPLETE  UNDERDETERMINED  SYSTEMS 

Realizing  from  what  is  demonstrated  above,  the  least  square  and/or  the  least  norm 
sense  in  which  Matlab  approaches  that  inverse  mathematical  problem,  it  is  essential  that  a 
thorough  understanding  of  how  this  linear  system  solution  is  depicted  to  the  dynamic 
system  in  terms  of  damage  detection,  is  possessed.  For  this  purpose,  an  FE  model  of  a 
cantilever  beam  consisting  of  ten  beam  elements  is  used.  First,  the  baseline  model  is  set 
and  the  modal  data  of  the  system,  i.e.,  natural  frequencies  and  mode  shapes,  are 
calculated,  along  with  the  stiffness  sensitivity  matrix  through  a  perturbation  of  one 
percent  reduction  to  the  rigidity  of  each  beam  element,  recursively.  Then,  another  FE 
model  is  built,  identical  to  the  baseline  with  the  only  difference  being  the  presence  of 
damage  somewhere  in  the  structure.  This  damage  is  represented  by  a  10  percent  rigidity 
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reduction  and  this  model  is  used  to  provide  simulated  modal  data  from  a  structure  that 
has  localized  damage.  We  will  refer  to  the  damage  as  the  "error,"  recognizing  that  these 
concepts  apply  equally  to  identification  of  model  error,  or  structural  damage.  The  main 
objective  is,  after  having  the  data  for  both  the  baseline  and  the  damaged  model  and  using 
the  stiffness  sensitivity  matrix,  to  examine  how  the  model  can  predict  the  true  element  in 
error,  which  is  the  design  variable  that  is  subject  to  the  10  percent  rigidity  reduction.  It  is 
noted  here  that  the  procedure  is  repeated  iteratively  for  all  the  beam  elements,  in  order  to 
assess  possible  errors  anywhere  in  the  model. 

1.  Using  Individual  Modes 

The  starting  point  of  this  analysis  is  the  use  of  only  one  mode  of  the  available 
data.  For  demonstration  purposes,  only  the  first  five  out  of  the  total  twenty  modes  will  be 
taken  into  consideration.  The  rows  of  the  stiffness  sensitivity  matrix  [ T ]  of  the  cantilever 

beam,  which  correspond  to  the  modes  of  the  system,  are  shown  in  the  Figures  1 1  and  12 
of  the  previous  chapter.  The  first  five  rows  are  reproduced  in  Figure  13.  Again  the  values 
are  normalized  to  the  maximum  value  of  each  row: 


33 


Mode  1 


Figure  13.  First  5  Rows  of  Stiffness  Sensitivity  Matrix 


The  numerical  representation  of  these  first  five  modes  is  given  below: 


0.0006 

0.0004 

0.0003 

0.0002 

0.0001 

0.0001 

0.0000 

0.0000 

0.0000 

0.0000 

0.0163 

0.0028 

0.0008 

0.0061 

0.0121 

0.0138 

0.0103 

0.0048 

0.0011 
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0.0904 

0.0094 

0.0691 

0.0778 

0.0170 

0.0221 

0.0977 

0.1083 

0.0407 

0.0025 

0.2518 

0.1633 

0.2888 

0.0390 

0.2834 

0.2793 

0.0435 

0.3602 

0.3188 

0.0291 

0.5564 

0.7459 

0.2616 

0.7508 

0.5051 

0.4994 

0.7337 

0.3113 

1.1042 

0.1708 

Now,  the  first  five  natural  frequencies  of  the  baseline  system  and  of  the  perturbed 
(damaged)  system,  and  the  first  five  entries  of  vector  jA&rj  formed  by  their  difference 
are  provided  from  the  FE  models. 

This  is  all  the  information  needed  for  solving  Equation  (2.31)  for  the  change  in 
the  design  variables  {ADF} .  The  following  figures  show  the  results  for  using  a  single 
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mode,  i.e.,  individual  modes  1  through  5,  respectively,  with  the  error  being  at  each 
element,  recursively.  In  all  the  figures,  the  y-axis  represents  the  percentage  error  and  the 
x-axis  the  number  of  the  element: 
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Figure  14.  Solution  Using  the  First  Mode  Only 
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Figure  15.  Solution  Using  the  Second  Mode  Only 
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Figure  16.  Solution  Using  the  Third  Mode  Only 
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Figure  17.  Solution  Using  the  Fourth  Mode  Only 
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Figure  18.  Solution  Using  the  Fifth  Mode  Only 

Keeping  in  mind  from  Figure  13  that  the  maximum  sensitivity  values  occur  at 
elements  1,  1,  8,  8,  9  for  modes  1,  2,  3,  4,  5,  respectively,  a  very  important  conclusion  is 
reached  for  the  FE  model  solution  using  only  individual  modes,  as  this  is  clearly  derived 

39 


from  the  study  of  Figures  14  through  18.  The  conclusion  here  is  that  no  matter  where  the 
actual  error  is,  the  calculated  result  depends  on  and  is  exclusively  led  by  the  highest 
sensitivity  value  of  the  individual  mode  used  for  the  solution.  This  covers  the  localization 
aspect  of  the  error  prediction.  Moreover,  it  is  observed  that  when  the  highest  sensitivity 
value  corresponds  to  the  element  being  actually  in  error,  then  not  only  will  the  location  of 
the  damage  be  detected,  but  also,  the  true  magnitude  of  the  error  will  be  revealed.  This  is 
the  reason  that  there  is  an  exact  correlation  of  the  calculated  error  to  the  actual  error  in 
cases  like:  first  element  in  error  using  first  or  second  mode,  eighth  element  in  error  using 
third  or  fourth  mode  and  ninth  element  in  error  using  fifth  mode.  It  can  be  predicted  at 
this  point  that  if,  for  example,  only  the  nineteenth  mode  is  used,  which  has  the  fifth 
element  with  the  highest  sensitivity  value,  only  when  the  fifth  element  actually  is  in  error 
will  there  be  a  positive  result,  which  can  be  verified  in  Figure  19. 
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Figure  19.  Solution  Using  the  Nineteenth  Mode  Only 
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2.  Using  More  Than  One  Mode 


The  next  step  of  the  current  analysis  will  be  the  examination  of  the  use  of  more 
than  one  mode,  starting  with  the  use  of  pairs  of  modes.  Consider  the  first  two  modes  and 
the  corresponding  stiffness  sensitivity  matrix  rows  from  the  two  top  plots  of  Figure  13. 
The  numerical  representation  of  them,  denoted  by  [l[  2  ]  and  using  four  significant  digits 
is: 


0.0006  0.0004  0.0003  0.0002  0.0001  0.0001  0.0000  0.0000  0.0000  0.0000 

0.0163  0.0028  0.0008  0.0061  0.0121  0.0138  0.0103  0.0048  0.0011  0.0000 


while  the  natural  frequencies  of  the  baseline  and  the  damaged  system  and  their  difference 
will  again  be  calculated  from  the  FE  models. 

Now,  solving  for  the  design  variables  using  the  first  two  modes,  we  get  the 
following  figure,  again  looping  over  all  the  elements  being  in  actual  error  recursively: 
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Figure  20.  Solution  Using  Modes  1  and  2 

As  shown  in  Figure  20,  the  solution  always  returns  a  vector  with  only  two  non¬ 
zero  values  (even  with  the  last  four  elements  being  in  error,  there  is  a  value  at  element  1 , 

but  it  is  relatively  minor  compared  to  that  of  the  sixth  element),  since  the  rank  of  the 
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sensitivity  matrix  used,  consisting  of  the  first  two  rows  of  the  full  sensitivity  matrix,  is  2. 
Furthermore,  these  non-zero  values  always  occur  at  elements  1  and  6  of  the  beam.  As 
explained  above,  it  is  no  surprise  that  a  positive  result  is  obtained  only  when  the  actual 
error  occurs  at  element  1  or  element  6  or  both,  as  will  be  shown  later.  In  order  to 
investigate  the  reason  for  this  trend,  the  procedure  for  the  solution  will  be  examined  in 
successive  steps  as  follows: 

Consider  again  the  stiffness  sensitivity  matrix  of  the  first  two  modes  shown 

above: 


0.0006  0.0004  0.0003  0.0002  0.0001  0.0001  0.0000  0.0000  0.0000  0.0000 

0.0163  0.0028  0.0008  0.0061  0.0121  0.0138  0.0103  0.0048  0.0011  0.0000 


This  is  a  case  of  an  underdetermined  problem,  which  as  discussed  previously  in 
this  chapter,  Matlab  solves  using  QR  decomposition  with  column  pivoting.  Column 
pivoting  means  that  the  columns’  norm  is  the  factor  taken  into  consideration  for  the 
pivoting  and  a  permutation  matrix  is  effectively  applied  to  the  columns  of  T.  It  is  of  vital 
importance  to  note  here  that  higher  modes  have  higher  sensitivities  and.  Therefore,  the 
higher  modes  tend  to  govern  the  nonn  of  the  sensitivity  matrix  columns.  Keeping  this  in 
mind  and  for  purposes  of  demonstration  it  can  be  considered,  without  any  loss  of 
generality,  that  a  partial  pivoting  is  adequate  to  explain  how  Matlab  reaches  this  solution 
for  this  linear  system. 


Moving  on  with  this  procedure,  the  largest  entry  that  is  considered  the  first  pivot 
is  entry  1  of  the  2nd  row. 


0.0006  0.0004 
0.0163  0.0028 


0.0003  0.0002 
0.0008  0.0061 


0.0001  0.0001 
0.0121  0.0138 


0.0000  0.0000 
0.0103  0.0048 


0.0000  0.0000 

0.0011  0.0000 


Performing  a  Gaussian  elimination  the  matrix  becomes: 


0  0.0003 
0.0163  0.0028 


0.0003  -0.0000 
0.0008  0.0061 


-0.0003  -0.0005 
0.0121  0.0138 


-0.0004  -0.0002 
0.0103  0.0048 


-0.0000  -0.0000 

0.0011  0.0000 
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The  next  pivot  will  be  chosen  from  the  first  row,  and  it  will  be  the  entry  with  the 
highest  absolute  value,  and  it  turns  out  to  be  entry  6: 


0  0.0003 
0.0163  0.0028 


0.0003  -0.0000 
0.0008  0.0061 


-0.0003  -0.0005 
0.0121  0.0138 


-0.0004  -0.0002 
0.0103  0.0048 


-0.0000  -0.0000 

0.0011  0.0000 


Performing  the  next  elimination,  the  matrix  becomes: 


0  0.0003 
0.0163  0.0130 


0.0003  -0.0000 
0.0091  0.0051 


-0.0003  -0.0005 
0.0019  0 


-0.0004  -0.0002 
-0.0006  -0.0004 


-0.0000  -0.0000 

-0.0001  -0.0000 


At  this  point,  it  can  arbitrarily  be  said  that  the  solution  is  governed  qualitatively 
by  the  highest  values  of  the  final  stage  of  the  elimination  of  the  rank-deficient  matrix,  as 
shown  above.  Comparing  Figure  20  to  this  final  stage,  in  terms  of  location,  the  solution 
for  the  design  variables  will  always  give  non-zero  values  at  the  pivots,  which  also  happen 
to  have  the  highest  values  in  each  row.  Again,  remember  that  this  is  not  absolutely  the 
way  Matlab  manipulates  its  solution  to  this  linear  problem,  but  it  serves  to  demonstrate 
the  main  ideas  of  this  result  and  to  illustrate  how  several  model  updating  methods  can  be 
properly  used. 


In  order  to  validate  what  is  stated  above,  modes  2  and  5  will  be  used  together  and 
an  attempt  to  predict  the  result  will  be  made.  The  distribution  of  these  particular 
sensitivity  rows  over  the  beam  is  shown  in  Figure  13,  while  their  arithmetic 
representation  in  a  common  matrix  is  given  below: 


0.0163  0.0028 
0.5564  0.7459 


0.0008  0.0061 
0.2616  0.7508 


0.0121  0.0138 

0.5051  0.4994 


0.0103  0.0048 
0.7337  0.3113 


0.0011  0.0000 
1.1042  0.1708 


The  first  pivot  will  be  entry  9  from  the  2nd  row: 


0.0163 

0.5564 


0.0028 

0.7459 


0.0008 

0.2616 


0.0061 

0.7508 


0.0121 

0.5051 


0.0138 

0.4994 


0.0103 

0.7337 


0.0048 

0.3113 


0.0011 

1.1042 


0.0000 

0.1708 
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From  this  point  already,  the  first  non-zero  solution  for  the  outcome  should  be  expected  at 
element  9.  The  elimination  based  on  that  pivot  gives: 


0.0158  0.0021  0.0006  0.0053  0.0116  0.0133  0.0096  0.0045  0  -0.0001 

0.5564  0.7459  0.2616  0.7508  0.5051  0.4994  0.7337  0.3113  1.1042  0.1708 


The  next  pivot,  that  is  the  1st  row  of  the  above  matrix,  will  be  entry  1 : 


0.0158  0.0021  0.0006  0.0053  0.0116  0.0133  0.0096  0.0045  0  -0.0001 

0.5564  0.7459  0.2616  0.7508  0.5051  0.4994  0.7337  0.3113  1.1042  0.1708 


Even  at  this  point,  the  second  non-zero  value  of  the  outcome  should  be  expected  at  the  1st 
element,  but  for  consistency  reasons,  the  next  elimination  is  presented: 


0.0158  0.0021 
0  0.6727 


0.0006  0.0053 
0.2422  0.5635 


0.0116  0.0133 
0.0965  0.0308 


0.0096  0.0045 
0.3955  0.1516 


0  -0.0001 
1.1042  0.1751 


So,  the  prediction  for  this  configuration  is  that  no  matter  where  the  actual  error  is, 
Matlab  will  return  a  result  with  non-zero  values  at  elements  1  and  9,  which  will  be 
positive  only  when  these  elements  are  actually  in  error,  either  one  or  both.  By  the  tenn 
positive  result,  the  correct  damage  detection  is  meant  in  a  qualitative  aspect  (error 
location),  which  might  sometimes  be  consistent  quantitatively  (error  magnitude).  Yet  the 
quantitative  aspect  cannot  be  proved,  although  in  the  cases  analyzed  in  this  subchapter,  it 
is  revealed. 

Solving  now  the  fundamental  Equation  (2.3 1)  for  the  design  variables  and  looping 
again  over  all  the  elements  that  are  in  error  recursively,  we  get: 
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Actual  Error 
Calculated  Error 


10%  Reduction  at  Elem  #  1 


10%  Reduction  at  Elem  #  2 


10%  Reduction  at  Elem  #  3 


10%  Reduction  at  Elem  #  4 


10%  Reduction  at  Elem  #  5 


10%  Reduction  at  Elem  #  8 


1  23456789  10 

10%  Reduction  at  Elem  #  9 


0. 

-0. 

0. 

-0. 


1  23456789  10 

10%  Reduction  at  Elem  #  10 


1  23456789  10 


Figure  2 1 .  Solution  Using  Modes  2  and  5 

It  is  easily  shown  in  Figure  21  that  what  is  stated  above  is  validated  exactly. 

In  order  to  validate  the  presented  procedure  in  a  more  generic  manner,  the  first 


five  modes  will  be  used  and  again  the  result  will  be  predicted  without  actually  solving 
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Equation  (2.31).  A  priori,  the  solution  expected  should  have  only  5  non-zero  values,  since 
the  rank  of  the  used  sensitivity  matrix,  shown  in  Figure  13  is  5.  In  this  case,  all  five  plots 
of  Figure  13  are  used  with  the  following  arithmetic  representation: 


'0.0006 

0.0004 

0.0003 

0.0002 

0.0001 

0.0001 

0.0000 

0.0000 

0.0000 

0.0000 

0.0163 

0.0028 

0.0008 

0.0061 

0.0121 

0.0138 

0.0103 

0.0048 

0.0011 

0.0000 

K:.]- 

0.0904 

0.0094 

0.0691 

0.0778 

0.0170 

0.0221 

0.0977 

0.1083 

0.0407 

0.0025 

0.2518 

0.1633 

0.2888 

0.0390 

0.2834 

0.2793 

0.0435 

0.3602 

0.3188 

0.0291 

0.5564 

0.7459 

0.2616 

0.7508 

0.5051 

0.4994 

0.7337 

0.3113 

1.1042 

0.1708 

The  first  pivot 

in  this  case  is  entry  9  of  the  5th  row: 

'0.0006 

0.0004 

0.0003 

0.0002 

0.0001 

0.0001 

0.0000 

0.0000 

0.0000 

0.0000 

0.0163 

0.0028 

0.0008 

0.0061 

0.0121 

0.0138 

0.0103 

0.0048 

0.0011 

0.0000 

fe]  = 

0.0904 

0.0094 

0.0691 

0.0778 

0.0170 

0.0221 

0.0977 

0.1083 

0.0407 

0.0025 

0.2518 

0.1633 

0.2888 

0.0390 

0.2834 

0.2793 

0.0435 

0.3602 

0.3188 

0.0291 

0.5564 

0.7459 

0.2616 

0.7508 

0.5051 

0.4994 

0.7337 

0.3113 

1.1042 

0.1708 

so  one  non-zero  value  should  be  expected  at  element  9.  After  the  elimination  based  on  the 
first  pivot,  the  matrix  becomes: 


'0.0006 

0.0004 

0.0003 

0.0002 

0.0001 

0.0001 

0.0000 

0.0000 

0  -0.0000 

0.0158 

0.0021 

0.0006 

0.0053 

0.0116 

0.0133 

0.0096 

0.0045 

0  -0.0001 

[?*]  = 

0.0699 

-0.0181 

0.0594 

0.0501 

-0.0017 

0.0037 

0.0707 

0.0968 

0  -0.0038 

0.0912 

-0.0521 

0.2133 

-0.1777 

0.1376 

0.1351 

-0.1683 

0.2704 

0  -0.0203 

0.5564 

0.7459 

0.2616 

0.7508 

0.5051 

0.4994 

0.7337 

0.3113 

1.1042  0.1708 

The  second  pivot  will  now  be  entry  8  of  the  4th 

row: 

'0.0006 

0.0004 

0.0003 

0.0002 

0.0001 

0.0001 

0.0000 

0.0000 

0  -0.0000 

0.0158 

0.0021 

0.0006 

0.0053 

0.0116 

0.0133 

0.0096 

0.0045 

0  -0.0001 

[?»]  = 

0.0699 

-0.0181 

0.0594 

0.0501 

-0.0017 

0.0037 

0.0707 

0.0968 

0  -0.0038 

0.0912 

-0.0521 

0.2133 

-0.1777 

0.1376 

0.1351 

-0.1683 

0.2704 

0  -0.0203 

0.5564 

0.7459 

0.2616 

0.7508 

0.5051 

0.4994 

0.7337 

0.3113 

1.1042  0.1708 

so  another  non-zero  value  should  be  expected  at  element  8.  After  the  elimination  based 
on  the  second  pivot,  the  matrix  becomes: 
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[4]= 


0.0006 

0.0004 

0.0003 

0.0002 

0.0001 

0.0001 

0.0000 

0 

0 

0.0000 

0.0142 

0.0029 

-0.0030 

0.0083 

0.0093 

0.0110 

0.0124 

0 

0 

0.0002 

0.0373 

0.0005 

-0.0169 

0.1137 

-0.0509 

-0.0447 

0.1309 

0 

0 

0.0034 

0.0912 

-0.0521 

0.2133 

-0.1777 

0.1376 

0.1351 

-0.1683 

0.2704 

0 

-0.0203 

0.5564 

0.7459 

0.2616 

0.7508 

0.5051 

0.4994 

0.7337 

0.3113 

1.1042 

0.1708 

It  must  be  said  here  that,  in  fact,  there  is  no  need  to  eliminate  entries  both  above 
and  below  the  pivots.  The  only  elimination  needed  for  the  purpose  of  this  demonstration 
is  of  the  remaining  rows  not  yet  containing  any  pivot,  which,  as  can  easily  be  shown, 
yields  to  the  form  of  an  upper  or  lower  triangular  matrix  after  the  appropriate  column 
shifting. 

The  third  pivot  will  be  entry  7  of  the  third  row: 


K»]= 


0.0006 

0.0004 

0.0003 

0.0002 

0.0001 

0.0001 

0.0000 

0 

0 

0.0000 

0.0142 

0.0029 

-0.0030 

0.0083 

0.0093 

0.0110 

0.0124 

0 

0 

0.0002 

0.0373 

0.0005 

-0.0169 

0.1137 

-0.0509 

-0.0447 

0.1309 

0 

0 

0.0034 

0.0912 

-0.0521 

0.2133 

-0.1777 

0.1376 

0.1351 

-0.1683 

0.2704 

0 

-0.0203 

0.5564 

0.7459 

0.2616 

0.7508 

0.5051 

0.4994 

0.7337 

0.3113 

1.1042 

0.1708 

making  element  7  likely  to  have  a  non-zero  value  at  the  result. 

Again,  after  the  elimination  of  the  remaining  rows  without  pivot,  the  matrix 
becomes: 


[*»]  = 


0.0006 

0.0004 

0.0003 

0.0002 

0.0001 

0.0001 

0 

0 

0  -0.0000 

0.0107 

0.0029 

-0.0014 

-0.0025 

0.0141 

0.0153 

0 

0 

0  -0.0001 

0.0373 

0.0005 

-0.0169 

0.1137 

-0.0509 

-0.0447 

0.1309 

0 

0  0.0034 

0.0912 

-0.0521 

0.2133 

-0.1777 

0.1376 

0.1351 

-0.1683 

0.2704 

0  -0.0203 

0.5564 

0.7459 

0.2616 

0.7508 

0.5051 

0.4994 

0.7337 

0.3113 

1.1042  0.1708 

and  now  the  next  pivot  is  entry  6  of  the  second  row,  which  makes  element  6  one  of  those 
with  a  non-zero  corresponding  value: 
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"0.0006 

0.0004  0.0003 

0.0002 

0.0001 

0.0001 

0 

0 

0  -0.0000" 

0.0107 

0.0029  -0.0014 

-0.0025 

0.0141 

0.0153 

0 

0 

0  -0.0001 

fo]  = 

0.0373 

0.0005  -0.0169 

0.1137 

-0.0509 

-0.0447 

0.1309 

0 

0  0.0034 

0.0912 

-0.0521  0.2133 

-0.1777 

0.1376 

0.1351 

-0.1683 

0.2704 

0  -0.0203 

0.5564 

0.7459  0.2616 

0.7508 

0.5051 

0.4994 

0.7337 

0.3113 

1.1042  0.1708 

After  the  last  elimination,  that  of  the  first  row,  the  matrix  becomes: 

"0.0005 

0.0004  0.0003 

0.0002 

0.0001 

0 

0 

0 

0  0.0000" 

0.0107 

0.0029  -0.0014 

-0.0025 

0.0141 

0.0153 

0 

0 

0  -0.0001 

[r*]  = 

0.0373 

0.0005  -0.0169 

0.1137 

-0.0509 

-0.0447 

0.1309 

0 

0  0.0034 

0.0912 

-0.0521  0.2133 

-0.1777 

0.1376 

0.1351 

-0.1683 

0.2704 

0  -0.0203 

0.5564 

0.7459  0.2616 

0.7508 

0.5051 

0.4994 

0.7337 

0.3113 

1.1042  0.1708 

In  this  last  matrix,  the  highest  value  of  the  first 

row  occurs  at  entry  1 : 

"0.0005 

0.0004  0.0003 

0.0002 

0.0001 

0 

0 

0 

0  0.0000" 

0.0107 

0.0029  -0.0014 

-0.0025 

0.0141 

0.0153 

0 

0 

0  -0.0001 

0.0373 

0.0005  -0.0169 

0.1137 

-0.0509 

-0.0447 

0.1309 

0 

0  0.0034 

0.0912 

-0.0521  0.2133 

-0.1777 

0.1376 

0.1351 

-0.1683 

0.2704 

0  -0.0203 

0.5564 

0.7459  0.2616 

0.7508 

0.5051 

0.4994 

0.7337 

0.3113 

1.1042  0.1708 

Having  reached  the  final  stage  of  this  procedure  using  the  first  five  inodes,  the 
predicted  solution  is  expected  to  have  non-zero  values  at  elements  1,  6,  7,  8,  9  and  these 
values  will  constitute  a  positive  result  when  the  actual  error  occurs  in  each,  or  any 
combination  of  these  elements. 

As  far  as  the  calculated  solution  from  the  solution  of  Equation  (2.31)  is 
concerned,  the  following  figure  is  generated,  which  again  validates  what  is  stated  above 
(again  where  the  existence  of  non-zero  values  at  elements  1,  6,  7,  8,  9  is  not  obvious  due 
to  the  very  low  value  relative  to  the  other  values): 
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Actual  Error 
Calculated  Error 


10%  Reduction  at  Elem  #  1 


Figure  22.  Solution  Using  Modes  1  through  5 

It  has  to  be  recalled  at  this  point  that,  although  the  quantitative  aspect  of  damage 
detection  (error  magnitude)  is  present  in  the  cases  examined  above,  it  is  not  proved  in  a 
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generic  way  for  all  situations.  It  will  later  be  shown  that  there  may  be  some  discrepancies 
which  distract  the  result  from  being  positive  as  far  as  the  error  magnitude  or  the 
exclusivity  of  the  location  are  concerned.  These  discrepancies  will  be  dealt  with  in  a  later 
chapter. 


3.  What  if  a  Fully  Determined  Square  Matrix  is  Available? 

From  what  is  stated  so  far,  a  very  important  conclusion  is  reached.  For  a  structure 
consisting  of  n  elements,  and  considering  one  design  variable  for  each  element,  if  m 
modes  are  used  for  the  sensitivity-based  analysis,  with  m<n,  then  the  maximum  number 
of  elements  that  can  be  positively  identified  to  be  in  error  is  k,  where  k  is  the  rank  of  the 
sensitivity  matrix,  which  at  most,  is  equal  to  m.  In  other  words,  if  it  is  assumed  that  there 
is  no  linear  dependency  between  the  rows  of  the  sensitivity  matrix,  then  the  maximum 
number  of  positive  damage-detected  elements  is  equal  to  the  number  of  the  modes  used 
for  the  analysis. 

Therefore,  it  can  be  concluded  that  if  damage  detection  throughout  the  structure  is 
desired  using  one  and  only  one  sensitivity  matrix,  then  this  matrix  must  not  only  be 
square,  but  also  full  rank.  This  is  the  condition  that  the  present  work  has  been  based  on.  It 
has  to  be  said  that  even  if  a  matrix  is  full  rank,  there  may  be  some  near-dependencies  that 
are  crucially  negative  for  the  desired  result.  This  is  why  the  enlargement  of  the  data  used 
must  first  be  pursued  and  then  the  optimal  configurations  contained  in  this  data  must  be 
found  and  used  for  globally  positive  structure  identification. 

The  following  figures  show  that  with  a  square  matrix,  comprised  of  either  the  first 
10  modes,  or  of  the  second  10  modes  of  the  structure  examined  previously  (cantilever 
beam),  potentially  good  results  are  obtained  (Figure  24)  because  all  the  elements  may 
have  a  non-zero  value  in  the  solution  vector.  Apparently,  this  is  only  necessary,  and  yet 
not  sufficient  condition  for  positive  damage  detection  (Figure  23),  since  different 
configurations  give  different  results  even  though  both  are  square  and  full  rank. 
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Figure  23.  Solution  Using  Modes  1  Through  10 
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Figure  24.  Solution  Using  Modes  1 1  Through  20 
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IV.  ENLARGEMENT  OF  TEST  DATA-  ENHANCEMENT  OF 
SENSITIVITY  MATRIX  DETERMIN ABILITY/  CONDITIONING 


In  the  previous  chapter,  the  need  for  the  construction  of  a  fully  determined 
sensitivity  matrix  for  the  detection  of  potential  damage  throughout  the  structure  was 
justified.  In  this  chapter,  several  methods  will  be  pursued  in  this  direction  until  the  point 
that  optimal  ABC  configurations  are  reached  and  prove  more  essential  than  other 
methods. 

A.  USING  EIGENVECTOR  SENSITIVITY  DATA 

In  the  theory  chapter,  a  method  was  presented  for  the  calculation  of  the 
eigenvector  sensitivity  matrix.  This  algorithm  is  applied  to  the  same  model  that  was 
analyzed  in  tenns  of  eigenvalue  sensitivity  and  the  results  are  provided  below. 

It  should  be  mentioned  at  this  point  that  even  if  data  from  eigenvector  are  used, 
resulting  in  an  improved  FEM  updating,  two  major  drawbacks  as  far  as  this  analysis  is 
concerned  cannot  be  ignored.  The  first  is  the  accuracy  of  the  test  measurements  needed  to 
accurately  identify  the  mode  shapes  of  a  structure  since  typical  measurement  errors  can 
be  caused  by  noise,  nonlinearities,  etc.  The  second  is  that  for  large  and  complex 
structures,  the  mode  shape  sensitivity  analysis  produces  a  large  quantity  of  data,  resulting 
in  a  more  expensive  problem  in  terms  of  calculation  time  and  memory  storage. 
Nevertheless,  a  possible  combination  of  the  two  aforementioned  discrepancies  may  cause 
the  whole  process  to  converge  less  rapidly  or  possibly  converge  into  incorrect  results  as 
described  in  [24]. 

For  this  analysis,  the  same  10  element  cantilever  beam  is  considered.  Using  the 

algorithm  described  in  Chapter  II,  a  full  eigenvector  sensitivity  matrix  is  generated.  The 

size  of  this  matrix  is  400x10.  The  number  of  columns  is  equal  to  the  number  of  the  design 

parameters,  which  again,  are  the  rigidities  corresponding  to  each  one  of  the  10  beam 

elements.  The  number  of  rows  corresponds  to  20  entries  related  to  each  DOF  of  each 

beam  element,  respectively,  which  entries  are  generated  for  each  one  of  the  20  modes  of 

the  system.  Obviously,  the  same  FE  model  that  was  used  for  the  derivation  of  the 
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eigenvalue  sensitivity  matrix  is  used  again.  The  distribution  of  the  eigenvector  sensitivity 
matrix  for  Mode  #1,  over  all  the  elements,  and  for  each  DOF  is  displayed  in  Figures  25 
and  26,  while  the  first  mode  of  the  eigenvalue  sensitivity  matrix  is  available  in  the  first 
plot  of  Figure  13.  Each  plot  is  normalized  to  the  maximum  value  of  each  row  for 
displaying  reasons.  It  should  be  mentioned  that  in  absolute  values,  the  magnitude  of  the 
eigenvector  sensitivities  is  much  less  than  that  of  the  eigenvalue  sensitivities,  as  shown  in 
the  arithmetic  representation  with  four  significant  digits  of  the  first  mode  of  each  one 
below,  and  this  may  play  a  significant  role  in  the  procedure  of  linear  system  solving: 

Teigyai, Model  =  10“3*[  0.6034  0.4388  0.3017  0.1928  0.1117  0.0565  0.0234  0.0070  0.0011  0.0000] 


T  =10-6* 

eigvec, Model 


-0.2217  0.0842 
-0.0721  0.0274 
-0.5550  0.0334 
-0.0399  -0.0423 
-0.7109  -0.1599 
-0.0129  -0.0228 
-0.7201  -0.2476 
0.0090  -0.0071 
-0.6131  -0.2520 
0.0259  0.0050 
-0.4191  -0.1946 
0.0380  0.0136 
-0.1653  -0.0948 
0.0459  0.0192 
0.1247  0.0306 
0.0502  0.0222 
0.4322  0.1682 
0.0519  0.0234 
0.7449  0.3094 
0.0522  0.0236 


0.0588  0.0381 
0.0191  0.0124 
0.2235  0.1447 
0.0353  0.0228 
0.2373  0.3081 
-0.0285  0.0313 
0.0992  0.3243 
-0.0179  -0.0237 
0.0172  0.1980 
-0.0098  -0.0187 
-0.0235  0.0971 
-0.0041  -0.0152 
-0.0362  0.0132 
-0.0004  -0.0130 
-0.0320  -0.0605 
0.0016  -0.0118 
-0.0198  -0.1295 
0.0023  -0.0113 
-0.0052  -0.1971 
0.0025  -0.0112 


0.0224  0.0115 
0.0073  0.0037 
0.0850  0.0436 
0.0134  0.0068 
0.1806  0.0925 
0.0183  0.0093 
0.3024  0.1545 
0.0221  0.0112 
0.2960  0.2263 
-0.0220  0.0126 
0.1702  0.1985 
-0.0201  -0.0198 
0.0536  0.0816 
-0.0189  -0.0192 
-0.0577  -0.0330 
-0.0183  -0.0190 
-0.1665  -0.1465 
-0.0180  -0.0189 
-0.2746  -0.2597 
-0.0180  -0.0189 


0.0048  0.0015 
0.0016  0.0005 
0.0183  0.0055 
0.0029  0.0009 
0.0387  0.0117 
0.0039  0.0012 
0.0646  0.0195 
0.0047  0.0014 
0.0944  0.0285 
0.0052  0.0016 
0.1269  0.0382 
0.0056  0.0017 
0.0915  0.0483 
-0.0156  0.0017 
-0.0016  0.0198 
-0.0155  -0.0098 
-0.0943  -0.0391 
-0.0154  -0.0098 
-0.1870  -0.0979 
-0.0154  -0.0098 


0.0002 

0.0000 

0.0001 

0.0000 

0.0009 

0.0000 

0.0001 

0.0000 

0.0019 

0.0001 

0.0002 

0.0000 

0.0032 

0.0001 

0.0002 

0.0000 

0.0046 

0.0002 

0.0003 

0.0000 

0.0062 

0.0002 

0.0003 

0.0000 

0.0078 

0.0003 

0.0003 

0.0000 

0.0095 

0.0003 

0.0003 

0.0000 

-0.0050 

0.0004 

-0.0042 

0.0000 

-0.0302 

-0.0025 

-0.0042 

-0.0007 
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Mode  1  DOF#1 
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Figure  25.  Mode  1  of  Eigenvector  Sensitivity  Matrix  for  DOF  1  Through  10 
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Mode  1  DOF  #11 


1  23456789  10 

Mode  1  DOF  #12 


1  23456789  10 

Mode  1  DOF  #14 


1  23456789  10 

Mode  1  DOF  #17 


Figure  26.  Mode  1  of  Eigenvector  Sensitivity  Matrix  for  DOF  1 1  Through  20 
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The  natural  frequency  and  the  eigenvector  differences  corresponding  to  the  first 
mode  can  be  derived  from  the  baseline  and  the  damaged  model  as  before.  In  this  case,  the 
^ Mode i  vector  will  have  a  dimension  of  20x1 . 

Two  approaches  will  now  be  made  to  the  identification  of  the  effectiveness  of  the 
use  of  eigenvector  sensitivity.  The  first  one  will  use  only  mode  shape  data  for  the 
analysis,  while  the  second  will  combine  the  same  modes  from  both  eigenvalue  and 
eigenvector  data.  It  is  quite  sufficient  to  use  only  the  first  mode  for  displaying  the  result 
for  both  cases  as  there  are  not  many  qualitative  variations  using  other  modes  as  well. 

Solving  the  fundamental  Equation  (2.31)  for  the  design  variables  using  only  the 
eigenvector  sensitivity  data  and  looping  over  all  the  elements  being  in  error  recursively 
yields  the  following  figure: 
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Figure  27.  Solution  Using  Mode  1  of  Full  Eigenvector  Sensitivity  Matrix  Only 

It  is  pretty  obvious  that  no  positive  result  can  be  extracted  from  this  approach, 
even  though  more  than  a  square  matrix,  an  over-determined  matrix,  is  used  for  the 
analysis. 
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However,  if  eigenvalue  sensitivity  data  is  used  along  with  mode  shape  data,  a 
totally  different  outcome  is  produced.  Using  only  the  first  mode,  as  seen  in  Figure  28,  this 
configuration  is  capable  of  predicting  damage  anywhere  in  the  structure: 
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Figure  28.  Solution  Using  Mode  1  of  Eigenvalue  and  Full  Eigenvector  Sensitivity  Matrix 
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The  same  positive  results  yield  if  only  the  translational  DOFs  are  considered  to  be 
measured  and,  thus,  be  considered  in  the  mode  shape  data,  as  is  shown  in  Figure  29. 
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Figure  29.  Solution  Using  Mode  1  of  Eigenvalue  And  All  Translational  DOF  of 

Eigenvector  Sensitivity  Matrix 
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In  spite  of  the  positive  results  shown  in  Figures  28  and  29,  it  is  realistically 
impossible  for  all  DOF,  or  even  all  the  translational  DOF  to  be  measured.  Therefore, 
spatially  incomplete  mode  shape  data  are  examined  in  the  following  two  ways.  The  first 
considers  DOFs  3,  7,  11,  15,  19  of  the  cantilever  beam  to  be  measured  and  provides 
accurate  mode  shape  data.  The  second  one  considers  DOFs  3,  5,  7,  9,  11,  13,  15,  17,  and 
19  to  be  measured.  Then,  the  first  mode  of  this  data  is  used  in  each  case  for  system 
identification  through  a  combined  eigenvalue  and  eigenvector  sensitivity  matrix,  which 
are,  respectively,  a  6x10  underdetermined  and  a  10x10  square  matrix.  The  solutions 
corresponding  to  these  cases  are  demonstrated  in  Figures  30  and  31.  The  figures  show 
that  these  are  not  configurations  that  give  positive  results  no  matter  where  the  true  error 
is.  Actually,  the  first  case  positively  identifies  the  error  in  six  out  of  ten  locations  (which 
in  terms  of  efficiency  is  the  best  it  could  do,  since  the  sensitivity  matrix  is  6x10,  as  was 
described  in  Chapter  III),  and  it  will  never  predict  the  error  in  the  remaining  four 
locations,  since  the  output  vector  has  only  six  fixed  non-zero  components.  On  the  other 
hand,  although  in  the  second  case  a  square  matrix  is  used,  it  is  not  able  to  predict  the 
correct  element  being  in  error  in  all  the  locations.  This  is  another  example  of  why  a 
square  matrix  is  a  necessary  (first  case,  Figure  30)  but  not  sufficient  (second  case,  Figure 
31)  condition  for  damage  detection  through  a  unique  sensitivity  matrix. 

Conclusively,  it  can  be  said  that  mode  shape  design  sensitivity  data  may  be  used 
for  the  augmentation  of  the  modal  frequency  data  used  for  the  sensitivity  based  FEM 
updating,  while  by  themselves  they  do  not  constitute  a  reliable  way  to  identify  errors  in 
the  structure.  Of  course,  considerable  thought  relative  to  the  right  choice  of  DOF  to  be 
measured  should  not  be  neglected,  because  different  configurations  have  quite  different 
effects  on  the  output  result  as  shown  in  the  cases  examined  above.  Beyond  that,  the 
drawbacks  referred  to  in  the  beginning  of  this  chapter  must  be  taken  into  account  in  terms 
of  the  availability  of  measurement  locations  and  of  the  accuracy  associated  with  those 
measurements. 
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Figure  30.  Solution  Using  Mode  1  of  Eigenvalue  and  Translational  DOFs  [3,  7,  11,  15, 
19]  of  Eigenvector  Sensitivity  Matrix  (Underdetennined) 


64 


Actual  Error 
Calculated  Error 


10%  Reduction  at  Elem  #  1 


Figure  31.  Solution  Using  Mode  1  of  Eigenvalue  and  Translational  DOFs  [3,  5,  7,  9,  11, 
13,  15,  17,  19]  of  Eigenvector  Sensitivity  Matrix  (Square) 
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B.  PARAMETER  SELECTION 

As  outlined  in  Chapter  I  under  the  background  section,  the  first  attempts  to 
override  discrepancies  that  are  generated  during  the  updating  procedure  were  focused  on 
what  is  called  parameter  selection.  This  method  can  be  viewed  in  different  variations, 
since  several  analyses  have  been  made  through  the  years  and  are  provided  in  the 
literature.  Some  of  them  will  briefly  be  described  in  the  following  paragraphs. 

1.  Subset  Selection 

The  parameter  subset  selection  method  chooses  a  subset  of  the  design  variables 
that  minimizes  a  penalty  function  based  on  the  norm  of  the  residual  of  Equation  (2.31)  as 
it  is  summarized  in  [5].  In  that  sense,  it  is  assumed  that  only  a  subset  of  the  parameters 
has  corresponding  independent  columns  of  the  sensitivity  matrix  and  the  optimum  subset 
is  chosen  through  an  iterative  procedure  to  produce  a  sub-optimal  solution,  commonly 
known  as  forward  selection  [3], 

The  first  step  is  to  find  a  single  column  among  the  columns  of  the  sensitivity 
matrix  that  best  represents  the  frequency  differences  vector  Am2.  Next,  the  combination 
of  two  columns  that  constitute  the  best  sub-basis  for  the  spanning  of  A  or  is  determined, 
while  the  first  parameter  selected  is  retained.  Then  a  Gram-Schmidt  orthogonalization  is 
performed  and  both  the  sensitivity  columns  and  the  Am2  are  replaced  by  the 
corresponding  vectors.  The  procedure  is  now  repeated  on  this  reduced  problem  for  the 
next  parameter  that  gives  the  smallest  residual.  An  iterative  process  is  then  produced  in 
order  to  combine  the  best  parameters  chosen  with  each  one  of  those  remaining  to  pursue 
the  best  data  representation.  It  should  be  noted  that  the  criterion  of  the  smallest  residual  is 
equivalent  to  the  smallest  angle  or  distance  of  the  vectors  or  the  vector  subspaces  existing 
or  created  respectively. 

The  subset  selection  problem  may  also  be  solved  by  using  an  exhaustive  search  of 
all  possible  subsets  that  minimize  the  penalty  function  and  eventually  choose  the  best 
subset  that  constitutes  the  global  minimum  of  all  the  cases.  However,  taking  into  account 
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that  the  number  of  parameters  and  also  the  size  of  the  subsets  may  be  large  for  reasonable 
realistic  problems,  the  computation  required  for  an  exhaustive  search  is  prohibitive  [5]. 

2.  Clustering  of  Design  Variables 

Another  approach,  similar  to  the  subset  selection  method,  is  that  dealing  with 
clustering  analysis.  In  this  approach,  the  columns  of  the  sensitivity  matrix  are  considered 
as  vectors  in  multi-dimensional  space.  For  example,  if  they  are  identical,  they  appear  as 
co-linear  vectors,  meaning  that  the  corresponding  parameters  affect  the  same  way  the 
modes  of  the  structure.  The  eigenvalue  sensitivity  matrix  for  a  large  FE  model  will 
contain  many  quasi-dependent  columns;  these  can  be  identified  via  the  clustering 
approach  and  may  be  combined  so  that  the  conditioning  of  the  problem  is  improved.  The 
most  similar  columns  are  added  together,  and  in  this  way,  the  information  available  from 
the  measurements  is  used  to  converge  the  model  by  the  smallest  adjustment  to  the 
updating  parameters  [9].  Again,  in  the  general  application  of  the  clustering  method,  a 
distance  criterion  is  needed  to  assess  the  similarity  of  different  objects  as  is  described  in 
[10]  and  then  as  is  also  described,  a  method  of  linking  the  similar  objects  within  the 
cluster  is  necessary.  Clustering  methods  have  been  successfully  demonstrated  in  [8],  [9], 
[10],  however,  it  has  to  be  said  that  if  during  the  updating,  the  sensitivities  columns 
change,  different  clusters  might  be  formed.  Nevertheless,  in  the  case  of  large  and 
complex  systems,  clustering  analysis  is  able  to  carry  out  the  task  of  efficiently  comparing 
many  possible  sensitivity  columns.  Together  with  sensitivity  analysis  and  subset 
selection,  it  provides  a  method  for  including  the  effects  of  changes  to  less  sensitive 
substructures  by  combining  those  columns  of  the  sensitivity  matrix  that  are  close  over  the 
bandwidth  of  the  modal  testing. 

3.  Grouped  Elements  Sensitivities 

A  different  manner  in  which  selecting  parameter  is  considered  is  presented  in  [7]. 
The  concept  of  total  sensitivity  is  introduced  as  the  absolute  sum  of  the  sensitivities.  The 
simple  observation  that  this  analysis  is  based  on  is  the  following.  Suppose  n  design 
variables  are  considered  and  they  are  about  to  be  merged  into  one  that  groups  in  a  way 
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the  effectiveness  of  the  selected  parameters  into  one  sensitivity  value.  Then,  if 
Ttot  represents  this  value  and  T{  represents  the  sensitivity  value  corresponding  to  each 
design  variable,  fori  =  1 :  n  ,  it  is  easy  to  see  that: 

Ttot-Tl+T2+-  +  Tn=fjTi  (4.1) 

1=1 

The  total  sensitivity,  as  defined  above,  will  be: 

\Ttot\  =  \Tl+T2+-  +  Tn\<fj\Ti\  (4.2) 

1=1 

The  equality  holds  only  when  the  signs  of  the  individual  sensitivity  values  of  the 
design  variables  have  the  same  sign.  This  leads  to  the  conclusion  that  by  grouping  design 
variables  into  larger  parameters,  the  number  of  parameters  can  be  reduced,  but  the  total 
sensitivity  decreases,  in  general.  From  this  basic  observation,  the  parameter  selection 
procedure,  sets  of  parameters  can  be  constructed  with  the  objective  that  their  effect  on  the 
residual  of  Equation  (2.31)  is  maximized,  meaning  that  the  residual  to  be  minimized  is 
more  sensitive  to  those  selected  parameters. 

This  procedure  is  accomplished  by  a  sequence  of  two  different  selection  phases 
[7].  The  first  phase  involves  merging  two  neighboring  elements  into  one  parameter,  as 
long  as  the  signs  of  their  sensitivity  values  are  the  same.  Then,  Equation  (4.2)  turns  out  to 
be  equality  and  identical  to  Equation  (4.1),  so  that  the  total  sensitivity  is  just  the  sum  of 
the  merged  parameters’  sensitivities.  This  phase  ends  when  no  neighboring  parameters 
have  the  same  sensitivity  sign. 

If  the  number  of  the  reduced  parameters  is  not  acceptable,  then  the  second  phase 
begins,  but  in  this  case,  the  sacrifice  of  total  sensitivity  is  unavoidable.  Obviously,  none 
of  the  neighboring  parameters  have  the  same  sign  of  sensitivities  and  by  merging  two 
neighboring  design  variables,  the  total  sensitivity  is  changed: 

ilfi+ld+l^bilfi+lr+f  (4.3) 

k= 1  k= 1 

k^i,j  k^i,j 

so  that  the  decrement  of  the  total  sensitivity  is  expressed  as: 
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(4.4) 


\Ti\  +  \Tj\-\Ti+Tj 

In  other  words,  the  sacrifice  for  reducing  one  parameter  in  terms  of  sensitivity,  in 
this  phase,  is  equal  to  Equation  (4.4).  This  phase  ends  when  some  criteria  that  have  been 
imposed,  such  as  the  number  of  design  variables  or  the  maximum  allowed  sacrifice,  or 
both,  are  met. 

In  this  work,  the  concept  of  grouping  neighboring  elements  is  used  to  enlarge  the 
available  data.  Since  a  cantilever  beam  is  examined,  and  as  seen  in  Figures  1 1  and  12,  the 
sensitivity  values  are  positive  for  all  the  design  variables  and  for  all  the  modes,  no  total 
sensitivity  sacrifice  is  required.  In  other  words,  the  first  phase  described  above  is 
sufficient  for  the  analysis.  So,  the  grouping  of  the  10  elements  of  the  cantilever  beam  into 
5  pairs  of  neighboring  elements  results  in  the  sensitivity  matrix  of  dimension  20x5  whose 
distribution  of  all  the  rows  is  presented  in  Figures  32  and  33.  Note,  that  in  the  following 
figures,  the  first  bar  corresponds  to  elements  1  and  2,  the  second  bar  corresponds  to  3  and 
4  and  so  on  until  the  fifth  bar,  which  corresponds  to  elements  9  and  10.  This  means  that, 
for  example,  the  highest  value  of  the  first  row  of  the  grouped  sensitivity  can  be  taken  to 
occur  at  either  the  first  or  the  second  element  of  a  10  element  cantilever  beam.  This  will 
prove  of  crucial  importance  later  on  when  the  concept  of  grouped  sensitivity  will  be  used 
to  enhance  the  use  of  ABC  for  this  analysis. 


69 


Mode  1 


1  2  3  4  5 

Mode  4 


1  2  3  4  5 

Mode  6 


1  2  3  4  5 

Mode  7 


Figure  32.  Grouped  Sensitivity  Modes  1  Through  10 
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Figure  33.  Grouped  Sensitivity  Modes  1 1  Through  20 
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c. 


THE  USE  OF  ABC 


Before  getting  into  the  details  of  using  ABC  as  an  efficient  way  of  enhancing  the 
sensitivity  matrix  determinability,  it  is  useful  to  mention  several  things  about  the  natural 
frequencies  of  the  ABC  systems  and  their  connection  to  the  baseline  system  natural 
frequencies,  with  a  short  reference  to  Cauchy’s  Interlace  Theorem. 

1.  Cauchy’s  Interlace  Theorem 

Interlace  theorem,  in  general,  deals  with  the  eigenvalues  of  a  principal  submatrix 
related  to  those  of  the  main  matrix.  According  to  the  analysis  made  in  [25],  assume  a 
matrix  A  of  dimension  nxn  and  the  submatrix  H  of  dimension  mxm  : 

with  A’s  eigenvalues:  ax<a-,  <a3---<an  and  H’s:  9X  <e2<e3- 

The  theorem  states  that  for  j  =  1, . . . ,  m  it  is: 

aj  ^  ^  aJ+n-m  and  equivalently,  a_u+n_m)  <  0_}  <  a.; . 

Alternatively,  for  k  -  1, 2, . . . ,  n  it  is: 

0k-n+m  ^  ak  ^  °i.  and  6  k  <  a_k  <  0_{k_n+m) . 

In  other  words,  6*  is  an  inner  bound  on  a±j  relative  to  A’s  spectrum.  For  example, 

[0  ^  0" 

and  A  =  e  0  1  Figure  34, 

J  [o  1  0 

below,  contains  a  schematic  representation  that  demonstrates  the  theorem  using  the 
values  of  the  above  example. 


consider  the  following  example,  where  H  = 


~H  B* 
B  H 
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Figure  34.  Schematic  Representation  of  Cauchy’s  Interlacing  Theorem 


The  actual  application  of  this  theorem  to  the  ABC  natural  frequencies  is  that  these 
frequencies  will  always  interlace  the  natural  frequencies  of  the  system  of  higher  order. 
This  means  that  if  a  baseline  system  is  considered,  like  the  cantilever  beam  of  the 
previous  analyses,  and  one  ABC  is  imposed  on  this  system,  then  the  eigenvalues  of  the 
ABC  system  will  interlace  the  eigenvalues  of  the  unconstrained  cantilever  beam.  So,  the 
first  eigenvalue  of  the  ABC  system  will  be  between  the  first  and  the  second  eigenvalue  of 
the  baseline  system  and  so  on.  If  a  second  ABC  is  imposed  on  the  system,  then  the  ABC 
system,  with  the  two  artificial  constraints,  will  have  natural  frequencies  that  will  interlace 
those  of  the  ABC  system  with  the  one  constraint.  Therefore,  in  terms  of  bandwidth 
limitations,  it  is  even  better  to  use  the  ith  mode  of  an  ABC  system  than  the  ith+1  mode  of 
the  baseline  system,  if  no  other  obvious  reason  for  the  use  of  each  is  present.  Figure  35 
may  be  useful  for  the  demonstration  of  this  interlacing: 

- ♦ - ♦ - ♦ - 


+ - ^ ^ ^ -  Eigenvalues  of  Baselin  e  with  1  ABC 


Eigenvalues  of  Baselin  e 


Ascending  Order 


Figure  35.  Schematic  Representation  of  Interlacing  Between  ABC  and  Baseline  Natural 

Frequencies 
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2.  Using  ABC 


Having  said  how  ABC  information  can  be  obtained  and  where  this  information  is 
relative  to  the  spectrum  of  the  baseline  system,  several  ABC  systems  are  used  to  facilitate 
the  error  prediction  in  the  cantilever  beam  that  is  examined  in  this  work.  All  these 
configurations  are  used  in  the  concept  of  building  an  at  least  square  sensitivity  matrix  that 
will  be  able  to  manipulate  data  from  the  FE  model,  and  from  the  modal  testing  of  the 
structure,  and  provide  a  positive  outcome  via  the  Matlab  linear  system  solution. 

The  manner  in  which  ABC  are  selected  is  based  on  the  principle  that  sensitivity  is 
increased  in  areas  near  the  artificially  imposed  constraint,  as  shown  in  [18].  This 
observation  can  help  since  the  sensitivity  distribution  throughout  the  structure  is  likely 
predictable,  at  least  for  the  fundamental  mode,  of  each  ABC  configuration.  Therefore,  the 
first  step  is  to  create  10  different  ABC  systems,  which  will  consist  of  the  baseline 
cantilever  beam  and  a  pin  imposed  on  each  one  of  the  10  unconstrained  translational 
DOF,  respectively.  All  the  data  needed  for  these  systems,  like  natural  frequencies,  mode 
shapes  and  sensitivities,  are  obtained  the  same  exact  way  as  for  the  baseline  system.  The 
first  mode  of  the  corresponding  sensitivities,  along  with  that  of  the  baseline  system,  is 
showed  in  Figure  36,  where  the  red  triangle  depicts  the  DOF  where  the  constraint  is 
imposed  and  all  the  values  are  nonnalized  to  the  maximum  value  of  the  first  mode  of 
each  system.  It  is  easily  noticeable  how  the  constraints  increase  the  sensitivity  near  the 
DOF  that  they  are  imposed  on,  as  long  as  these  DOF  are  not  near  the  free  end  of  the 
cantilever.  As  one  can  see,,  the  three  last  configurations  tend  to  simulate  the  second  mode 
of  the  baseline  system,  which  is  not  strange  considering  the  values  of  their  natural 
frequencies  (close  to  the  second  baseline  eigenvalue)  and  their  mode  shapes  (qualitatively 
close  to  the  second  baseline  mode  shape). 
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Figure  36.  First  Mode  of  Sensitivity  Matrix  of  Baseline  (Plot  1)  and  of  Each  ABC 

System  (Plots  2-11) 


Next,  the  solution  to  the  fundamental  Equation  (2.31)  is  pursued  using  two 


approaches.  The  first  one  uses  all  the  configurations  described  in  Figure  36,  thus 
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constructing  an  overdetermined  11x10  sensitivity  matrix,  while  the  second  uses  only  the 
ABC  configurations  (Plots  2-11  of  Figure  36),  thus  constructing  a  square  matrix.  These 
solutions  are  demonstrated  in  Figures  37  and  38,  where  the  prediction  of  this  model  for 
all  the  elements  being  in  error  recursively  is  provided. 
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Figure  37.  Solution  Using  Mode  1  of  Baseline  System  and  of  Each  ABC  System 
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Figure  38.  Solution  Using  Mode  1  of  Each  ABC  System  Only 

The  two  solutions  provided  in  the  above  figures  are  quite  similar  since  the 
sensitivity  matrices  used  in  each  case  do  not  have  many  differences  in  terms  of  sensitivity 
distribution  along  the  structure,  and  since  the  frequencies  difference  vector  used  is  almost 
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the  same.  What  is  also  obvious,  is  that  the  use  of  these  ABC  configurations  is  quite 
encouraging  since,  at  no  cost,  a  relatively  good  solution  is  reached  with  the  actual  error 
being  predicted  in  all  cases,  both  in  location  and  in  magnitude,  the  only  drawback  being 
the  presence  of  bias  considerably  comparable  to  the  true  error  in  some  of  the  cases  (when 
actual  element  in  error  is  1,4,  5,  6,  7,  8).  Towards  the  direction  of  eliminating  this  bias, 
other  configurations  will  be  pursued  later  in  this  thesis. 

A  very  important  observation  that  is  made  from  studying  Figure  36  and  imagining 
these  plots  as  the  rows  of  a  sensitivity  matrix  is  that,  as  was  described  before,  the  last 
three  rows  seem  almost  identical.  This  means  that  they  add  a  degree  of  linear  dependency 
between  them,  which  does  not  provide  the  linear  system  with  any  useful  additional 
information,  but  also  constitutes  a  discrepancy  in  the  whole  procedure  of  the  linear 
system  solution  as  was  stated  in  the  introductory  chapter.  Therefore,  only  one  of  these 
three  configurations  may  essentially  be  used  as  a  means  of  adding  data  to  the  system 
while  the  other  two  must  be  replaced  by  other  more  useful  configurations.  It  can  be  seen 
that  even  by  excluding  the  configuration  that  corresponds  to  a  constraint  imposed  on  the 
eighth  translational  DOF  (plot  9  of  Figure  36),  which  is  the  one  closer  to  the  second 
mode  of  the  baseline  system,  and  thus  using  a  square  sensitivity  matrix,  the  solution  of 
Figure  39  is  achieved.  This  solution  is  somehow  similar  to  those  of  Figures  37  and  38, 
but  yet  not  best  since  there  is  still  a  degree  of  linear  dependency  between  the  last  two 
rows. 
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Figure  39.  Solution  Using  Mode  1  of  Baseline  and  of  Each  ABC  System  Excluding  the 
One  That  Constrains  the  Eigth  Translational  DOF 

Moreover,  it  can  be  observed  that  up  to  the  point  that  the  similar  rows  described 


above  are  excluded,  the  remaining  sensitivity  matrix  presents  a  trend  that  reveals  a  very 
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large  degree  of  linear  independency  and  well-conditioning  and  which  is  the  milestone  of 
the  present  work.  This  is  the  concept  of  a  diagonally  dominant  or  a  tridiagonal  dominant 
matrix  that  is  described  in  the  following  subchapter  and  which  will  be  sought  through 
baseline,  ABC  and  grouped  element  sensitivities  configurations.  This  kind  of  matrix  is  a 
kind  of  banded  matrix  in  tenns  of  having  its  larger  values  on  the  main  diagonal  and  on 
the  first  diagonal  below  and/or  above  the  main  diagonal.  In  sensitivity  terms,  it  means 
that  there  exists  a  sensitivity  concentration  around  all  the  different  elements  of  the 
structure,  uniquely  identified  in  each  row  of  the  sensitivity  matrix.  This  kind  of  matrix 
will  be  called  diagonally  dominant  banded  for  the  purpose  of  this  thesis  since,  to  the 
author’s  knowledge;  there  is  no  official  mathematical  definition  that  can  describe  it. 

The  first  attempt  at  constructing  such  a  sensitivity  matrix  using  modes  from  the 
baseline  system  and  from  the  ABC  systems  described  so  far  is  made  using  first’  fifth  and 
twentieth  modes  of  the  baseline  system  (plots  1,  5  of  Figure  1 1  and  plot  10  of  Figure  12) 
and  the  first  mode  of  each  ABC  system  that  imposes  a  constraint  to  translational  DOF  1 
through  7  respectively  (plots  2  through  8  of  Figure  36).  That  way,  an  effort  to  achieve 
high  sensitivity  values  to  all  the  elements  of  the  cantilever  beam  in  different  sensitivity 
rows  is  made,  along  with  keeping  out  rows  that  do  not  provide  any  useful  information  for 
the  linear  system  solution,  such  as  the  three  last  plots  of  Figure  36.  The  sensitivity  matrix 
then  built  is  demonstrated  in  Figure  40  and  the  solution,  with  respect  to  the  design 
variables,  is  provided  in  Figure  41.  A  second  attempt  uses  the  same  configurations 
described  in  the  first,  replacing  the  fifth  mode  of  the  baseline  with  the  second  mode  of  the 
ABC  system  that  constrains  the  eighth  translational  DOF.  The  sensitivity  matrix 
constructed  and  the  corresponding  solutions  are  provided  in  Figures  42  and  43, 
respectively. 
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Figure  40.  Sensitivity  Matrix  of  Mode  1,  5,  20  of  Baseline  (Plots  l,  9,  10)  and  Mode  1  of 
Each  ABC  System  That  Constrains  Translational  DOF  1  Through  7  (Plots  2-8) 
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Figure  41. 
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Solution  Using  Mode  1,  5,  20  of  Baseline  and  Mode  1  of  Each  ABC  System 
That  Constrains  Translational  DOF  1  Through  7 
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Figure  42.  Sensitivity  Matrix  of  Mode  l,  20  of  Baseline  (Plots  1,10)  And  Mode  1  of 

Each  ABC  System  That  Constrains  Translational  DOF  1  through  7  (Plots  2-8)  and 
Mode2  of  ABC  System  That  Constrains  Translational  DOF  8  (Plot  9) 
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Figure  43.  Solution  Using  Mode  1,  20  of  Baseline  and  Mode  1  of  Each  ABC  System 
That  Constrains  Translational  DOF  1  Through  7  and  Mode2  of  ABC  System  That 

Constrains  Translational  DOF  8 


84 


As  seen  in  the  previous  figures,  the  solutions  provided  by  using  these 
configurations  are  quite  accurate,  both  in  location  and  in  magnitude,  as  far  as  the  error 
prediction  is  concerned.  Of  course,  these  examples  deal  with  a  very  simple  structure,  for 
which  several  ABC  configurations  with  not  only  one  constraint  can  be  used,  as  long  as 
natural  frequencies  of  those  systems  are  within  the  frequency  spectrum  available. 
Speaking  of  the  available  frequency  spectrum,  it  should  also  be  noted  that  modes  higher 
than  those  implemented  by  the  nature  of  each  problem,  which  in  general  are  a  few  lower 
modes,  cannot  be  used  as  it  is  done  with  the  twentieth  mode  of  the  baseline  system  in  the 
above  examples.  Knowing  this,  and  that  the  use  of  ABC  systems  may  not  be  unlimited  as 
mentioned,  it  is  shown  later  in  this  work  how  other  methods  of  enhancing  the  test  data, 
like  grouped  element  sensitivities,  can  be  combined  with  the  use  of  ABC  systems,  aiming 
at  the  construction  of  a  diagonally  dominant  banded  sensitivity  matrix. 

3.  Diagonally  Dominant  Matrix 

A  short  description  of  a  diagonally  dominant  matrix  is  provided  in  [26]. 
According  to  that  reference,  a  matrix  is  said  to  be  diagonally  dominant  if  in  every  row  of 
the  matrix,  the  magnitude  of  the  diagonal  entry  in  that  row  is  larger  than  or  equal  to  the 
sum  of  the  magnitudes  of  all  the  other  (non-diagonal)  entries  in  that  row.  The  concept  can 
be  extended  to  column  entries.  More  precisely,  the  matrix  A  is  diagonally  dominant  if: 

\au\>  ^  |a.;  for  all  i ,  where  atj  denotes  the  entry  in  the  ith  row  and  jth  column. 

j*i 

It  is  very  important  to  note  here  that  if  a  linear  system  is  examined  where 
elementary  row  operations,  such  as  switching  rows,  are  allowed  either  through  a 
pennutation  matrix  or  just  by  switching  the  corresponding  entries  of  the  right  hand  side 
vector,  instead  of  examining  the  criterion  of  diagonal  dominance  for  the  diagonal  entry, 
the  maximum  value  of  each  row  can  be  considered.  Then,  if  the  matrix  under 
consideration  has  a  maximum  value  at  different  entries  for  each  row  covering  its  whole 
length,  with  suitable  row  switching,  it  can  take  the  form  of  a  diagonally  dominant  matrix. 
Implementing  this  concept  in  the  global  sensitivity  matrix  used  for  model  updating,  if 
appropriate  configurations  for  a  structure  can  be  found  and  combined  together  such  that 
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in  each  row,  different  elements  are  subject  to  the  highest  sensitivity  values  while  the  rest 
of  the  values  satisfy  the  criterion  described  above,  this  matrix  can  be  said  to  be  a 
diagonally  dominant  one.  Furthermore,  if  the  strict  mathematical  definition  of  diagonal 
dominance  is  not  satisfied,  but  there  is  a  significant  concentration  of  high  sensitivity 
values  around  the  diagonal  entries,  then  this  matrix  can  be  called  a  diagonally  dominant 
banded  matrix,  as  stated  previously  in  this  chapter. 

The  main  issue  regarding  this  type  of  matrix  is  that  its  properties  make  it  of 
particular  use  in  several  cases.  It  has  been  proved  that  convergence  is  guaranteed  if  a 
diagonally  dominant  matrix  is  used  in  many  methods,  like  Jacobi  or  Gauss  Seidel,  for 
solving  linear  systems.  Also,  no  partial  pivoting  is  necessary  for  a  strictly  column 
diagonally  dominant  matrix  when  perfonning  Gaussian  elimination.  By  the  Gershgorin 
circle  theorem,  a  strictly  (or  irreducibly)  diagonally  dominant  matrix  is  non-singular.  This 
result  is  known  as  the  Levy-Desplanques  theorem  [26].  These  applications  are  of  major 
importance  in  FE  model  updating  since  they  guarantee  that  there  is  going  to  exist  a 
solution  to  the  linear  system,  which  can  be  solved  well-conditioned  and  with  no  concern 
about  dependencies  that  may  diverge  the  solution.  The  only  concern  is  finding 
configurations  that  can  sufficiently  add  data  to  the  underdetermined  problem,  aiming  at 
creating  a  diagonally  dominant  sensitivity  matrix. 

D.  MANIPULATING  GROUPED  ELEMENT  SENSITIVITIES 

It  has  been  mentioned  above  that  grouping  neighboring  elements  by  simply 
adding  their  sensitivities  values,  with  no  sacrifice  of  the  total  sensitivity  as  long  as  they 
have  the  same  sign,  is  an  effective  way  of  selecting  parameters.  Figures  32  and  33  show 
how  the  grouped  element  sensitivities  are  constructed  for  each  mode  of  the  cantilever 
beam  examined,  using  pairs  of  two  neighboring  elements.  The  next  thing  to  be 
considered  is  how  these  configurations  can  essentially  be  used  to  enlarge  the  available 
data  and  to  assist  in  the  effort  of  constructing  a  diagonally  dominant  banded  matrix. 
Doing  so  can  lead  to  a  sensitivity  matrix,  which  combines  baseline,  ABC  and  grouped 
element  sensitivities  information. 
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Consequently,  it  must  be  shown  that  grouped  element  sensitivities  data  can 
replace  other  data  in  the  global  sensitivity  matrix.  At  first  sight,  this  seems  to  be  a 
difficult  task,  since  from  what  is  shown  in  Figures  32  and  33,  if  only  one  miscellaneous 
mode  is  examined,  its  length  will  be  five,  corresponding  to  the  five  pairs  formed  from  the 
grouped  elements.  But  the  global  sensitivity  matrix  where  this  grouped  element 
sensitivity  mode  is  to  be  used  has  length  10,  corresponding  to  the  number  of  individual 
elements  of  the  structures  which  are  not  grouped.  Note  here  that  the  desirable  use  of 
additional  data  is  towards  the  enhancement  of  the  detenninability  of  the  sensitivity 
matrix,  and  not  towards  the  design  variables  increase,  and  that  can  be  performed  by  using 
grouped  element  sensitivities  as  extra  rows  for  the  global  sensitivity  matrix.  Therefore, 
there  must  be  a  way  found  for  this  five-component  row  to  be  pennuted  to  a  10  component 
row  and  then  added  to  the  global  sensitivity  matrix,  without  any  loss  of  its  properties. 

For  that  purpose,  the  twentieth  mode  will  be  considered  in  for  the  concept  used 
for  that  pennutation  to  be  explained.  The  highest  value  of  this  mode  occurs  in  the  last 
pair  of  the  cantilever  beam  elements,  which  are  elements  9  and  10.  Assume  that  this 
value  is  denoted  by  T^9andW.  Then,  the  rate  of  change  of  the  sum  of  the  rigidities  of 
elements  9  and  10,  with  respect  to  the  change  of  the  twentieth  natural  frequency  will  be 
equal  toT™9andl0 .  In  other  words: 


y»  20 

gr' 9  and  10 


_  rri  20  .  rji  20  _ 


A DV 


gr. 9  and  10 


ADF,„  +ADV, 


(4.5) 


From  Equation  (4.5)  it  can  be  seen  that  this  particular  sensitivity  value  can  be 
used  to  characterize  and  identify  the  potential  change  in  the  design  variable  of  element  9 
or  10  or  both,  due  to  the  change  of  the  twentieth  natural  frequency  of  the  system.  That 
leads  to  the  conclusion  that  T™9andl0  can  be  used  as  the  sensitivity  value  corresponding  to 

one  of  the  elements  9  or  10,  while  the  other  will  have  zero  sensitivity  value,  when  both 
elements  are  considered  in  a  sensitivity  row.  This  can  be  equivalently  expressed  by  the 
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fact  that  if,  for  example,  one  of  the  two  rigidities  of  elements  9  and  10  is  reduced  by  10 
percent,  that  is  either  ADVeigor  A DVelw  is  0.1,  then  the  total  decrease  to  their  group  will 

be  10  percent  too. 

This  conclusion  allows  a  wide  range  of  possible  manipulations  of  the  grouped 
element  sensitivities  data  in  order  to  be  appropriately  used  as  enhancement  measures  of 
the  global  sensitivity  matrix  determinability.  Figure  44  demonstrates  two  possible  ways 
of  how  the  twentieth  mode  of  the  grouped  element  sensitivity  can  be  used  in  the  global 
sensitivity  matrix  (of  length  10)  for  providing  a  sensitivity  row  with  highest  value  at  both 
element  9  and  element  10,  by  simply  adding  zeros  to  the  elements  that  are  not  required  to 
have  the  sensitivity  value  desired. 

Permuted  Grouped  Element  Sensitivity  of  Mode  20  with  Highest  Value  at  Element  9 


Permuted  Grouped  Element  Sensitivity  of  Mode  20  with  Highest  Value  at  Element  10 


Figure  44.  Mode  20  of  Grouped  Element  Sensitivity  Permuted  to  Attribute  Highest  Value 
Either  to  Element  9  (Plot  1)  or  to  Element  10  (Plot2) 
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E. 


COMBINING  BASELINE,  ABC  AND  GROUPED  ELEMENT 
SENSITIVITIES  FOR  A  DIAGONALLY  DOMINANT  BANDED  MATRIX 


The  ultimate  purpose  of  this  work  is  to  combine  baseline,  ABC  and  grouped 
element  configurations  in  a  an  essential  manner  into  a  global  sensitivity  matrix,  which 
will  be  able  to  predict  change  to  the  design  variables  of  a  system,  given  the  difference 
between  the  natural  frequencies  of  the  system  to  measure  (potentially  damaged)  and  those 
of  the  original  (undamaged)  model.  Having  said  that,  the  configurations  of  Figure  40  will 
be  used  again,  that  is  modes  1,5,20  of  the  baseline  system  and  mode  1  of  each  ABC 
system  that  constrains  the  translational  DOF  1  through  7,  with  two  different  alterations. 
In  the  first,  mode  5  of  baseline  system  will  be  replaced  by  the  twentieth  mode  of  grouped 
element  sensitivities  permuted  to  have  the  highest  value  at  the  ninth  element  (plot  1  of 
Figure  44),  while  in  the  second  mode  20  of  baseline  system  will  be  replaced  by  the 
twentieth  mode  of  grouped  element  sensitivities  permuted  now  to  have  the  highest  value 
at  the  tenth  element  (plot  2  of  Figure  44).  The  concept  of  these  replacements  is  to  prove 
that  modes  of  the  baseline  system  that  contribute  to  the  diagonally  dominant  banded 
matrix  with  high  sensitivity  values  at  certain  elements,  can  be  mutually  replaced  by 
modes  of  grouped  element  sensitivities  that  have  high  values  at  the  same  elements, 
respectively.  It  is  even  better,  under  the  principles  instituted  for  the  diagonally  dominant 
banded  matrix,  if  the  replacing  modes,  apart  from  only  high  value  at  certain  elements, 
have  a  significant  sensitivity  concentration  around  the  element  with  the  highest 
sensitivity  value,  but  this  is  something  that  is  achievable  depending  on  the  nature  of  each 
problem. 

Considering  the  first  case  mentioned  above,  the  configurations  used  are  presented 
in  terms  of  sensitivity  distribution  per  row  in  Figure  45,  and  the  solution  reached  using 
these  configurations  is  depicted  in  Figure  46. 
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Figure  45.  Sensitivity  Matrix  of  Mode  1, 20  of  Baseline  (Plots  1,10)  and  Mode  1  of  Each 
ABC  System  That  Constrains  Translational  DOF  1  Through  7  (Plots  2-8)  and  Mode  20 
of  Grouped  Element  Sensitivities  With  Highest  Value  at  Element  9  (Plot  9) 
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Figure  46.  Solution  Using  Mode  1,  20  of  Baseline  and  Mode  1  of  Each  ABC  System 
That  Constrains  Translational  DOF  1  Through  7  and  Mode  20  of  Grouped  Element 
Sensitivities  With  Highest  Value  At  Element  9 
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It  is  obvious  from  Figure  46  that  this  combination  of  configurations  works  fine  for 
predicting  the  error  when  it  occurs  in  elements  1  through  8,  but  it  is  not  positive  for  errors 
near  the  free  end  of  the  cantilever  beam.  It  has  to  be  mentioned  here  that  even  if 
sensitivity  matrix  is  full  rank,  there  is  a  discrepancy  regarding  the  frequency  difference 

vector  |  Aft)2 1 .  This  discrepancy  is  related  to  the  use  of  the  same  values  in  its  last  two 

entries,  since  the  twentieth  natural  frequencies  of  the  baseline  and  the  experimental 
model  are  the  same  for  both  the  original  and  the  grouped  element  model.  This  is  a  very 
important  observation  because  it  leads  to  the  conclusion  that  there  is  a  limitation  imposed 
when  grouped  element  sensitivities  configurations  are  used  along  with  baseline  ones. 
This  limitation  defines  that  the  same  modes  of  baseline  and  grouped  element  sensitivities 
cannot  be  used  together  in  a  global  sensitivity  matrix,  because  the  subject  of  using 
repeated  eigenvalues  is  raised  and  this  produces  discrepancies  to  the  linear  system 
solution.  This  limitation  can  be  further  applied  to  the  use  of  ABC  systems  that  have 
natural  frequencies  very  close  to  those  of  the  baseline  system,  as  was  discussed  for  the 
first  eigenvalue  of  the  ABC  system  that  constrains  the  eighth  translational  DOF  (plot  9  of 
Figure  36)  in  comparison  with  the  second  eigenvalue  of  the  baseline  system. 

Considering,  now  the  second  case  where  the  twentieth  mode  of  the  baseline 
system  is  replaced  by  the  twentieth  mode  of  the  grouped  element  sensitivities  permuted 
to  have  highest  sensitivity  value  at  the  tenth  element,  Figure  47  is  generated  by  the 
sensitivity  distribution  along  the  rows  of  the  sensitivity  matrix  created  by  these 
configurations.  No  repeated  eigenvalues  limitation  is  present  and  Figure  48  is  produced 
to  demonstrate  the  solution  given  from  the  linear  system  using  this  matrix. 
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Figure  47.  Sensitivity  Matrix  of  Mode  1,  5  of  Baseline  (Plots  1,  9)  and  Mode  1  of  Each 
ABC  System  That  Constrains  Translational  DOF  1  through  7  (Plots  2-8)  and  Mode  20  of 
Grouped  Element  Sensitivities  With  Flighest  Value  at  Element  10  (Plot  10) 
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Figure  48.  Solution  Using  Mode  1,  5  of  Baseline  and  Mode  1  of  Each  ABC  System  That 
Constrains  Translational  DOF  1  Through  7  and  Mode  20  of  Grouped  Element 
Sensitivities  With  Highest  Value  at  Element  10 
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Figure  48  provides  a  great  result  as  far  the  location  and  the  magnitude  of  the  error 
is  concerned,  wherever  in  the  structure  it  may  occur.  It  clearly  demonstrates  that  baseline, 
ABC  and  grouped  element  sensitivities  configurations  can  be  essentially  combined  into 
one  global  sensitivity  matrix  for  the  structure,  which  has  the  properties  of  a  diagonally 
dominant  banded  matrix  as  described  previously  in  this  work.  With  proper  attention  to  all 
the  limitations  mentioned  before,  it  is  possible  to  reveal  any  potential  damage  in  a 
structure  or  to  update  an  FE  model. 

In  addition  to  what  is  concluded  above,  the  same  concept  can  be  extended  to  cases 
where  actual  error  occurs  in  more  than  one  element  of  the  structure.  For  this  purpose, 
Figure  49  is  generated  and  includes  the  solution  using  the  configurations  depicted  in 
Figures  41,  47  and  45,  respectively,  with  actual  error  occurring  simultaneously  in 
elements  1,  3  and  10.  Again,  it  is  proved  that  optimal  configurations  are  those  who  adopt 
the  concept  of  a  diagonally  dominant  sensitivity  matrix,  but  also  considering  issues 
generated  by  using  similar  modes  of  different  systems. 
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Error  using  modes  1 ,5,20  of  baseline  and  1  st  of  each  2-8  ABC  syztem 


Error  using  modes  1 ,5  of  baseline, 1st  of  each  2-8  ABC  syztem  and  20th  of  grouped  at  10th  element 


Error  using  modes  1 ,20  of  baseline, 1  st  of  each  2-8  ABC  syztem  and  20th  of  grouped  at  9th  element 
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Figure  49.  Solution  For  Actual  Error  in  Elements  1,3,  10  Using:  Plot  1 -Configurations  of 
Figure  41,  Plot  2-Configurations  of  Figure  47,  Plot  3 -Configurations  of  Figure  45 
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V.  CONCLUSIONS/RECOMMENDATIONS 


CONCLUSIONS 

1 .  The  natural  frequencies  of  the  ABC  systems  can  easily  be  identified  from 
a  square  experimentally  measured  FRF  matrix.  These  natural  frequencies 
will  interlace  the  natural  frequencies  of  the  baseline  system  in  accordance 
with  the  Cauchy  interlace  theorem.  The  eigenvalue/eigenvector 
sensitivities  of  these  ABC  systems  can  be  calculated  using  the  same  FE 
model  as  that  used  to  generate  the  baseline  system  sensitivities. 

2.  There  is  an  exact  and  direct  correlation  between  the  stiffness  sensitivity 
matrix  and  the  flexural  properties  of  the  structure.  This  means  that  strain 
energy  modes  have  the  same  distribution  in  the  structure  as  their 
corresponding  stiffness  sensitivity  modes. 

3.  For  overdetennined  full  rank  problems,  the  use  of  Matlab  operators  “\” 
and  “pinv”  both  return  the  same  least  square  solution  to  the  linear  system 
created  by  the  fundamental  Equation  (2.31),  which  minimizes  the 
Euclidean  norm  of  the  residual,  once  with  QR  decomposition  and  once 
with  the  SVD  method. 

4.  For  underdetermined  or  rank  deficient  problems,  the  solutions  obtained 
from  the  use  of  those  two  operators  are  not  the  same.  In  this  case,  “pinv” 
gives  a  least  nonn  solution,  in  the  sense  that  the  length  of  the  solution 
vector  is  minimum,  while  “\”  returns  a  basic  solution  with  the  significance 
of  having  a  number  of  non-zero  values  equal  to  or  less  than  the  rank  of  the 
sensitivity  matrix. 

5.  Using  only  one  mode  for  the  identification  of  the  structure,  the  “\” 
operator  gives  a  solution  that  has  one  non-zero  value  corresponding  to  the 
element  with  the  highest  sensitivity  value  ,  while  using  the  “pinv”  operator 
returns  a  solution  that  is  scaled  by  the  sensitivity  distribution  of  the  mode 
used  and  eventually  by  the  strain  energy  of  the  same  mode. 

6.  The  use  of  “\”  turns  out  to  be  more  effective  for  system  identification 
because  it  provides  a  more  transparent  result  in  tenns  of  error  localization 
and  it  makes  clear  that  the  use  of  modes  with  locally  concentrated 
sensitivity  values  will  reveal  errors  at  the  corresponding  elements  where 
errors  exist. 

7.  Solving  underdetennined  systems,  that  is  using  fewer  modes  than  the 
number  of  elements  potentially  in  error,  will  identify  potential  damage  in, 
at  most,  as  many  elements  as  the  number  of  modes  used.  This  means  that 
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an  error  can  be  uniquely  identified  anywhere  in  a  structure  only  if  a  square 
full  rank  sensitivity  matrix  is  used  that  is  of  a  size  equal  to  the  number  of 
elements  potentially  in  error. 

8.  Mode  shape  design  sensitivity  data  may  sometimes  be  used  for  the 
augmentation  of  the  modal  frequency  data  used  for  the  sensitivity  based 
FEM  updating,  as  long  as  the  availability  and  suitability  of  the 
measurement  locations  is  appropriate,  and  the  accuracy  associated  with 
those  measurements  is  within  the  requirements  tolerance. 

9.  A  useful  means  of  compensating  for  an  insufficient  amount  of  test  data  is 
the  grouped  element  sensitivities  permuted  in  such  a  way  as  to  match  the 
sensitivity  matrix  length  by  simply  adding  zeros  in  elements  that  are  not 
required  to  have  a  sensitivity  value,  and  the  total  sensitivity  of  the  group 
being  assigned  to  the  remaining  member  of  that  group.  Grouped  element 
sensitivities  are  not  useful  when  the  same  modes  of  the  baseline  system 
are  used,  because  issues  of  repeated  eigenvalues  and  perhaps  of  mode 
quasi-dependencies  are  raised. 

10.  Baseline,  ABC  and  grouped  element  sensitivities  configurations  can  be 
used  together  in  one  sensitivity  matrix  with  the  only  limitation  being  the 
use  of  modes  from  different  systems  that  might  be  close  to  each  other, 
hence  creating  quasi-dependencies  and  close-to-repeated  eigenvalues 
issues. 

11.  If  baseline,  ABC  and  grouped  element  sensitivities  are  used  to  create  a 
“global”  square  sensitivity  matrix,  which  clearly  concentrates  high 
sensitivity  values  around  every  element  in  the  structure  in  different  modes, 
under  the  concept  of,  as  described,  a  “diagonally  dominant  banded” 
matrix,  then  error  can  be  identified  anywhere  in  this  structure  both 
quantitatively,  in  tenns  of  location,  and  qualitatively,  in  terms  of 
magnitude. 

RECOMMENDATIONS 

1 .  Verily  that  the  concept  of  the  “diagonally  dominant  banded”  matrix  can  be 
applied  to  the  identification  of  errors  associated  with  mass/density 
parameters. 

2.  Create  an  algorithm/program  that  chooses  the  best  combination  from  a 
variety  of  ABC  and/or  grouped  element  sensitivities  sets,  along  with 
baseline  sensitivities,  into  a  sensitivity  matrix  that  can  identify  a  parameter 
error  regardless  of  its  location.  This  should  be  approached  under  the 
constraint  of  the  availability  of  modes  in  a  certain  test  bandwidth  and  the 
criterion  of  minimum  computational  time  relative  to  the  size  and  the 
complexity  of  the  structure.  If  the  conditions  set  above  make  this 
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procedure  prohibitive,  find  combinations  of  configurations  that  create 
underdetermined  sensitivity  matrices,  capable  of  damage  prediction,  each 
one  in  different  areas/elements  of  the  structure,  as  some  regions  of  a 
structure  may  be  more  neuralgic  than  others. 
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